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ABSTRACT 


The  stability  characteristics  of  various  compact  fourth-  and  sixth-o^der  spatial  opera¬ 
tors  are  assessed  using  the  theory  of  Gustafsson,  Kreiss  and  Sundstrom  (G-K-S)  for  the 
semi-discrete  Initial  Boundary  Value  Problem  (IBVP).  These  results  are  then  generalized  to 
the  fully  discrete  case  using  a  recently  developed  theory  of  Kreiss.  In  all  cases,  favorable 
comparisons  are  obtained  between  G-K-S  theory,  eigenvalue  determination,  and  numerical 
simulation.  The  conventional  definition  of  stability  is  then  sharpened  to  include  only  those 
spatial  discretizations  that  are  asymptotically  stable  (bounded,  Left  Half-Plane  eigenval¬ 
ues).  It  is  shown  that  many  of  the  higher-order  schemes  which  are  G-K-S  stable  are  not 
asymptotically  stable.  A  series  of  compact  fourth-  and  sixth-order  schemes,  which  are  both 
asymptotically  and  G-K-S  stable  for  the  scalar  case,  are  then  developed. 
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1.  INTRODUCTION 


Recently,  higher-order  numerical  methods  have  seen  increasing  use  in  the  Direct  Numeri¬ 
cal  Simulation  (DNS)  of  the  Navier-Stokes  equations.  Although  they  do  not  have  the  spatial 
resolution  of  Spectral  methods,  they  offer  significant  increases  in  accuracy  over  conventional 
second-order  methods.  They  can  be  used  on  any  smooth  grid,  and  do  not  have  an  overly 
restrictive  CFL  dependence  as  compared  with  the  0(N~2)  CFL  dependence  observed  in 
some  Spectral  methods  on  finite  domains.  In  addition,  they  are  generally  more  robust  and 
less  costly  than  Spectral  methods.  The  issue,  of  the  relative  cost  of  higher-order  schemes 
(accuracy  weighted  against  physical  and  numerical  cost)  is  a  very  complex  issue,  ultimately 
depending  on  what  features  of  the  solution  are  sought  and  how  accurately  they  must  be 
resolved.  In  any  event,  the  further  development  of  the  underlying  stability  theory  of  these 
schemes  is  important. 

The  state  of  1  igher-order  temporal  discretizations  is  well  developed  and  relies  greatly 
on  the  existing  Ordinary  Differential  Equation  (ODE)  literature.  Higher-order  spatial  dis¬ 
cretizations  are  wed  documented  in  the  literature,  as  well.  For  example,  entire  classes  of 
centered  explicit  spatial  schemes  are  described  in  the  text,  of  Kopal  [1].  In  practice,  compact 
schemes  (methods  where  both  the  solution  and  its  derivatives  are  treated  as  unknowns  and 
are  solved  for  simultaneously)  are  more  accurate  than  optimal  explicit  schemes  (nth  order 
schemes  involving  n+1  grid-points),  and  have  gained  much  attention  recently  for  use  in  DNS. 
The  fundamental  ideas  of  compact  schemes  as  well  as  derivation  techniques  can  be  found  in 
the  work  of  Vichnevetsky  [2],  and  will  not  be  pursued  in  this  work.  The  primary  difficulty 
in  using  higher-order  schemes  is  finding  stable  boundary  schemes  that  preserve  their  formal 
accuracy.  It  is  well  known  that  for  a  hyperbolic  system  to  preserve  formal  spatial  accuracy, 
an  (N)th  order  inner-scheme  must  be  closed  with  at  least  an  (N  —  \)lh  order  boundary  scheme 
[3].  To  date,  many  higher-order  inner-schemes  are  used  with  lower-order  boundary  schemes, 
because  no  stable  higher-order  formulations  are  known.  The  formal  accuracy  of  these  for¬ 
mulations  is  thus  reduced  to  one  order  more  than  the  Boundary  Condition  (BC)  accuracy, 
and  it  is  questionable  whether  the  additional  work  incurred  in  the  higher-order  inner-scheme 
is  justified. 

Determining  the  numerical  stability  of  a  fully  discrete  approximation  for  a  linear  hy¬ 
perbolic  partial  differential  equation  is  a  difficult  task.  For  the  “Cauchy  Problem"  on  an 
infinite  domain  (  —  oc.oo),  standard  techniques  based  on  Fourier  methods  generally  provide 

neces^ry  conditions  for  stability  of  the  numerical  scheme.  For  the  Initial-Boundary- 
Valuc  Problem  (IBVP)  on  the  semi-infinite  domain  [0,  oo),  or  the  finite  domain  [-1,1],  Fourier 
techniques  are  not  straight  forward  to  apply  and  do  not  provide  sufficient  conditions  for  nu¬ 
merical  stability.  To  address  these  issues,  Osher  [4]  and  Kreiss  [5],  and  later  Gustafsson  et 


al.  [6],  developed  stability  analysis  techniques  based  on  normal  mode  analysis.  Their  work 
(generally  referred  to  as  G-K-S  stability  theory)  established  conditions  that  the  inner  and 
boundary  schemes  must  satisfy,  to  ensure  stability.  The  G-K-S  theory  states  that  hyperbolic 
systems  are  assured  stability  if  no  eigenvalues  or  generalized  eigenvalues  exist  for  the  IBVP. 
Trefethen  [7]  further  clarified  the  physical  meaning  of  the  G-K-S  condition  by  noting  that 
the  concept  of  stability  at  a  boundary  could  be  related  to  the  group  velocity  of  the  boundary 
scheme,  specifically,  whether  it  is  carrying  energy  into  or  out  of  the  numerical  domain. 

One  of  the  weak  points  of  fully  discrete  G-K-S  theory  has  been  the  great  complexity 
m  applying  it  to  higher  md<_i  uum^Kcai  schemes.  Raising  the  spatial  or  temporal  accuracy 
generally  increases  the  complexity  of  the  stability  polynomials  (the  order  increases,  giving 
rise  to  more  roots)  which  govern  the  stability  of  the  numerical  approximation.  In  multi-stage 
time  discretization  schemes  (e.g.  Runge-Kutta  schemes  with  three  or  more  stages),  where 
boundary  conditions  must  be  applied  at  intermediate  levels,  the  stability  polynomials  that 
must  be  tested  at  each  boundary  arc  nearly  insurmountable  using  analytic  techniques.  One 
can  greatly  simplify  the  analysis  by  addressing  the  semi-discrete  problem  as  a  method-of-lines 
IBVP,  rather  than  the  fully  discrete  problem.  The  underlying  G-K-S  theory  for  the  semi- 
discrete  problem  was  initially  developed  by  Strikwerda  [8],  He  showed  that  by  discretizing 
space  and  leaving  time  continuous,  the  necessary  and  sufficient  conditions  for  method-of-lines 
IBVP  stability  are  analogous  to  those  governing  the  stability  of  the  fully  discrete  case. 

The  precise  connection  between  the  semi-discrete  stability  bounds  and  those  obtained  in 
the  fully  discrete  analysis  is  not  always  straight  forward.  Recently,  however,  Kreiss  ct  al.  [9] 
have  shown  that  under  very  weak  conditions,  stability  of  the  semi-discrete  approximation 
infers  stability  in  the  fully  discrete  approximation  if  specific  Runge-Kutta  time  marching 
schemes  are  used.  Therefore,  one  can  rely  on  the  semi-discrete  G-K-S  theory  to  assess  the 
stability  with  R-K  integration,  of  various  higher-order  spatial  discretization  operators,  thus 
simplifying  the  calculations  considerably. 

The  emphasis  of  this  work  is  to  apply  semi-discrete  G-K-S  theory  to  several  higher- 
order  spatial  operators  for  the  IBVP.  Since  compact  methods  naturally  lend  themselves 
to  fewer  implementational  difficulties  at  the  boundaries,  they  will  be  the  primary  focus. 
Stable  boundary  formulations  which  preserve  the  formal  accuracy  of  the  inner-scheme  will 
be  presented  for  spatial  derivative  operators  of  up  to  sixth-order. 

2.  MODEL  EQUATION  FOR  IBVP 

Under  the  assumption  of  locally  frozen  coefficients,  the  equations  governing  conservation 
of  mass,  momentum  and  energy  in  one  spatial  dimension  can  be  transformed  into  a  system 
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(1) 


of  hyperbolic  equations,  having  the  form: 


dU_ 

dt 


+  F;  x  >  0,  t  >  0 


where  “A”  is  a  diagonal  matrix  with  real  eigenvalues,  F  is  a  source  term,  and  boundary  and 
initial  data  are  supplied  in  the  form 


U,(0,t)  +  TUI\0,t)=g(t),  U(x,0)  =  f(x) 


(2) 


where  T  is  a  matrix  describing  the  boundary  conditions,  and  A  has  been  divided  into  its 
right-  and  left-going  characteristics  to  determine  the  boundary  conditions.  The  problem  is 
said  to  be  well-posed  if  the  solution  U(x,t)  depends  smoothly  on  the  initial  and  boundary 
data.  Our  focus  in  this  work  wull  be  on  the  scalar  form  of  eqn  (1-2),  where  the  matrix  “A” 
is  a  negative  real  constant  “a”  and  with  the  source  term  F  =  0.  This  simplification  can 
be  justified  since  stability  of  a  numerical  scheme  on  the  scalar  equation  implies  stability  of 
the  system  if  boundary  conditions  are  imposed  in  a  characteristic  form  (ex.  see  Gottlieb  et 
al.  [10]).  A  further  necessary  modification  is  that  the  semi-infinite  spatial  domain  must  be 
truncated  to  a  finite  domain.  For  simplicity,  we  shall  assume  that  the  physical  domain  is 
confined  to  the  interval  0  <  x  <  1  for  all  times. 


3.  SEMI-DISCRETIZATION 

In  this  work,  the  numerical  discretization  of  cqn  (1)  will  be  accomplished  by  two  separate 
and  independent  steps.  The  spatial  derivatives  will  first  be  approximated  with  appropriate 
formulas,  leaving  what  is  generally  referred  to  as  a  semi-discretization.  The  numerical  solu¬ 
tion  will  then  be  advanced  in  time  in  a  Method-of- Lines  approach,  using  a  stable  temporal 
scheme.  We  begin  by  dividing  the  continuous  domain  [0,1]  into  N  uniform  intervals  of  width 
Ac  where  N Ax  =  1.  The  continuous  derivative  Ux  is  then  replaced  with  a  finite-difference 
representation  involving  the  functional  values  Uj  at  the  discrete  points.  A  system  of  ODE’s 
results  having  the  form 

Al  =  A/+^  j  =0,  ...,N  (3) 

where 

M+Vj  =  £  «  mkjVj+k  (4) 

k=-L, 

and  “L”  and  “R”  are  the  width  of  the  stencil  extending  to  the  left  and  right  of  grid-point 
j,  respectively.  Note  that  rn^Lj  and  are  functions  of  “j”  since  there  is  no  reason  to 


assume  that  the  same  stencil  will  be  used  at  each  grid-point.  For  example,  consider  the 
case  where  the  same  spatial  stencil  is  used  at  every  interior  grid-point  in  the  domain.  We, 
therefore,  choose  a  particular  =  L  and  Rj  =  R.  (Obviously,  the  choice  of  L,  R,  and 
nij  should  coincide  with  a  scheme  that  is  stable  for  the  Cauchy  problem  on  the  infinite 
domain).  This  scheme  can  only  be  used  for  L  <  j  <  N  —  R  without  the  stencil  protruding 
through  the  boundary.  Thus,  exactly  L  -f  R  additional  formulas  have  to  be  defined  near  the 
boundaries.  Since  there  is  only  one  physical  boundary  condition,  L  +  R  -  1  of  the  schemes 
are  strictly  numerically  motivated.  These  schemes  are  generally  referred  to  as  numerical 
boundary  schemes  (NBS). 

Noting  that  the  physical  boundary  condition  g(t)  will  bi.  imposed  at  the  grid-point  j— 0, 
eqn  (3)  can  be  rewritten  as 

=  MVj  +  Bj9{t);  j  =  l ,  ...,iV  (5) 

where  M  is  an  N  x  N  matrix,  and  Bj  is  a  vector  of  dimension  N,  describing  the  dependence 
of  the  "jth"  scheme  on  the  boundary  data.  The  matrix  M  usually  is  diagonal  of  order  L  +  R 
+  1  for  most  explicit  methods,  but  can  in  general  be  full.  With  little  loss  of  generality,  g(t) 
can  be  chosen  to  be  zero.  The  exact  solution  of  the  semi-discretization  described  by  eqn  (5) 
for  homogeneous  boundary  data  becomes 

Vj(t)  =  f(xj)expMt;  j  =  1,  (6) 


Note  that  all  the  boundary  information  is  incorporated  directly  into  the  matrix  M.  and  that 
the  stability  of  the  numerical  scheme  depends  directly  on  the  properties  of  the  matrix,  not 
just  on  what  spatial  discretization  was  chosen  for  the  inner-scheme. 

It  is  instructive  to  clarify  these  points  involving  boundary  closures  with  examples  of  an 
explicit  and  an  implicit  spatial  scheme.  We  choose  to  concentrate  on  schemes  having  at 
least  fourth-order  spatial  accuracy  since  most  of  the  difficulties  associated  with  high-order 
stencils  are  not  observed  with  second-order  schemes.  The  first  is  an  explicit  5-point  scheme 
reported  by  Gary  [11],  and  later  shown  by  Strikwerda  [8]  to  be  G  K-S  stable.  The  scheme  is 
uniformly  fourth-order  in  space.  The  spatial  discretization  is  accomplished  with  the  stencils 
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The  terms  in  eqn  (3)  take  the  form 
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where  M+  —  j£^M+.  Note  that  four  NBS’s  are  required  to  close  the  numerical  scheme,  of 
wrhich  only  one  will  be  overwritten  by  a  physical  boundary  condition.  Accounting  for  the 
boundary  condition  at  the  grid-point  j=0  in  eqn  (8),  results  in  expressions  for  eqn  (5)  of  the 
form 
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where  M  — 


M  and  D  = 


-D.  It  is  evident  that  the  matrix  M  is  the  N  x  N  submatrix 


12A*  ‘  ^  —  12AxJ 

of  the  matrix  M+,  corresponding  to  1  <  i,j  <  N. 

The  implicit  fourth-order  example  is  that  of  a  compact  discretization  in  which  the  nu- 
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merical  approximation  to  is  made  in  the  following  form 
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Noting  the  structure  of  the  spatial  operator  and  the  fact  that 
eqn  (3)  is  often  rearranged  into  the  form 


dV, 

dt 


av, 

1  dx 


M2  Vj ;  j  =  0,  ...,N 


,  for  this  scheme, 


(11) 
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where  A/+  =  (A/*)  1 A// •  The  resulting  matrix  expression  becomes 
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where  A/2+  =  ^A/2.  Imposing  the  boundary  information  at  grid- point  j=0  reduces  the 
matrix  by  one  order  and  yields  matrix  equations  of  the  form 
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where  /i2  =  2 .  Note  that  because  a  matrix  multiplication  was  involved  in  determining 

the  spatial  operator  for  grid-points  1  <  j  <  /V,  matrix  M  is  not  a  simple  submatrix  of  matrix 
Af+  as  it  was  with  the  explicit  scheme.  The  spatial  operator  described  by  eqn  (13)  is  referred 
to  as  being  compact  because  both  Mi  and  A/2  only  involve  three  grid-points  (except  near 
the  boundary  in  A/2).  Bccauc..  of  this  compactness,  L  =  R  =  1  ,  and  only  two  NBS  must 
be  defined  (one  of  which  includes,  but  is  not  replaced  by,  the  physical  boundary  condition). 
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Expressing  this  algorithm  in  the  form  of  eqn  (5)  yields 

M  =  B  =  M{XB2.  (14) 

It  is  apparent  that  while  the  algorithm  is  implementationally  compact,  the  resulting,  M  is  a 
full  N  x  N  matrix. 

Working  with  the  scalar  form  of  eqn  (1),  with  the  boundary  conditions  posed  at  x=0, 
has  allowed  us  to  simplify  the  semi-discretization  defined  by  eqn  (3)  to  that  of  eqn  (5).  Note 
that  if  the  governing  equation  would  have  been  Ut  4-  aUx  =  F  with  boundary  data  at  x 
=  1,  the  two  previously  mentioned  schemes  could  have  easily  accommodated  this  change 
and  still  produced  a  stable  algorithm.  For  a  system  of  hyperbolic  equations  with  variable 
coefficients,  one  does  not  know  a  priori  the  sign  of  the  eigenvalues  of  the  matrix  A  in  eqn  (1). 
Therefore,  the  solution  is  advanced  in  time  for  j=0,N,  followed  by  a  characteristic  decom¬ 
position  in  which  the  physical  boundary  conditions  are  imposed  at  either  j=0  or  j=N  such 
that  the  resulting  system  is  well-posed.  This  requirement  imposes  a  constraint  on  the  class 
of  allowable  matrices  il7+  which  can  be  conveniently  used  to  obtain  stable  discretizations. 
Only  central-difference  schemes  have  this  desirable  feature;  that  of  being  stable  for  a  >  0,  or 
a  <  0.  Assuming  that  a  central-difference  operator  is  used  throughout  the  interior  domain, 
then  eqn  (4)  requires  that  R  =  L  NBS’s  be  defined  at  each  boundary.  The  structure  of  these 
NBS’s  must  be  such  that  the  resulting  scheme  be  stable  for  either  the  inflow  (a  <  0;  j=l,N) 
or  the  outflow  (a  >  0;  j=0,N-l).  This  can  only  be  accomplished  if  the  same  methodology  is 
used  to  derive  the  NBS’s  at  each  end  of  the  domain.  Therefore,  the  NBS’s  at  each  respective 
boundary,  are  asymmetric  mirror  images  of  each  other.  The  resulting  matrix  M+  has  the 
following  property 

M+  =  —PM+P  (15) 

where  P  is  the  permutation  operator  defined  by  p,,j  =  0, 0  <  i,j  <  N  ,  for  i  ^  N  —  and 
Pi,j  =  1, 0  <  i,j  <  N  ,  for  i  =  N  -  j.  Note  that  PP  =  I,  the  identity  matrix.  The  matrix  M 
is,  therefore,  the  N  x  N  submatrix  of  M+. 

We  now  have  a  well  defined  class  of  spatial  operators  which  are  acceptable  for  our  pur¬ 
poses.  They  are  high-order  central  difference  schemes  with  boundary  implementations  which 
are  imposed  asymmetrically.  For  future  reference,  we  shall  define  a  convenient  nomencla¬ 
ture  to  describe  these  matrices.  The  matrix  M+  shall  be  described  as  (NBSi,''  NBSl  - 
CD  -  NBSr ,  •  •  • NBSs ),  where  “CD”  is  the  order  of  the  central  difference  operator  used 
in  the  interior,  and  NBSj  is  the  order  of  the  NBS  used  to  close  the  scheme  at  the  points 
next  to  the  boundary.  For  example,  the  explicit  uniformly  fourth-order  scheme  represented 
by  matrix  eqn  (8)  is  denoted  by  the  nomenclature  (4, 4-4-4, 4),  where  the  “-4-”  denotes  the 
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inner-scheme  being  approximated  with  a  fourth-order  stencil,  and  the  symmetric  “4,4*  de¬ 
notes  fourth-order  stencils  at  j=0,l,  and  j=N-l,N.  The  three-point  compact  scheme  described 
by  eqn  (12)  is  denoted  (4-4-4).  Again,  the  “-4-”  refers  to  the  inner-scheme  accuracy,  and 
the  symmetric  “4”s  indicate  closure  on  the  boundaries  of  fourth-order  accuracy.  There  is 
no  ambiguity  in  the  nomenclature  between  the  compact  and  explicit  schemes,  since  only  the 
compact  scheme  can  retain  fourth-order  inner  accuracy  with  one  NBS  at  each  boundary. 

4.  STABILITY  OF  THE  IBVP 

The  eigenvalues  of  matrix  M  from  eqn  (5),  resulting  from  higher-order  finite-difference 
approximations  to  Ux,  tend  to  align  along  the  imaginary  axis  in  complex  conjugate  pairs. 
To  time  advance  eqn  (5)  efficiently,  the  time  discretization  algorithm  should  include  in  its 
stability  domain  a  large  portion  of  the  imaginary  axis.  Conventional  Runge-Kutta  (R-K) 
time  advancement  algorithms  of  third-  or  fourth-order  are  well  suited  for  semi-discretizations 
of  hyperbolic  equations,  and  are  the  only  method  considered  in  this  study.  In  particular, 
the  standard  fourth-order  method  of  Kutta  (ex.  see  Gear  [12])  will  be  used  in  this  study 
because  of  the  1)  fourth-order  non-linear  accuracy,  2)  the  large  stability  envelope,  and  3)  the 
low  storage  requirements. 

It  is  desired  to  know  whether  a  given  higher-order  spatial  discretization  is  stable  for  this 
time  advancement  scheme.  Using  conventional  G-K-S  theory  for  the  fully  discrete  IBVP  for 
fourth-order  R-K  time,  and  second-order  space  discretization  would  involve  polynomials  of 
eighth-order  in  k  to  solve  at  the  boundaries.  Closed  form  solutions  would  be  difficult  to 
obtain  under  these  circumstances,  and  would  get  more  complicated  with  increasing  spatial 
accuracy.  The  stability  analysis  can  be  greatly  simplified  by  relying  on  three  fundamental 
theorems  of  stability  analysis  which  are  valid  under  the  conditions  proposed  in  this  study. 
We  will  discuss  each  briefly,  but  for  further  clarification  suggest  consulting  the  original  works. 
The  essential  elements  of  the  theorems  forming  the  basis  of  this  work  are: 

Theorem  1:  G-K-S  theory  (fully  discrete  [6]  or  semi-discrete  [8])  asserts  that  to  show 
stability  for  the  finite  domain  problem,  it  is  sufficient  to  show  that  the  inner-scheme  is 
Cauchy  stable  on  (— oo,  oo)  ,  and  that  each  of  the  two  quarter-plane  problems  is  stable  using 
normal  mode  analysis.  Thus,  the  stability  of  the  finite-domain  problem  is  broken  into  the 
summation  of  three  simpler  problems. 

Theorem  2:  For  each  quarter  plane  problem  arising  in  theorem  1,  a  necessary  and 
sufficient  condition  for  stability  of  the  initial  boundary  value  problem  is  that  there  exist  no 
eigensolution.  This  theorem  is  true  for  either  the  fully- discrete  case  [6],  or  the  semi-discrete 
case  [8]. 
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The  algebraic  complexity  involved  in  showing  stability  of  the  IBVP  is  dramatically  re¬ 
duced  in  the  semi-discrete  case,  since  time  remains  continuous.  Ultimately,  numerical  stabil¬ 
ity  is  a  fully-discrete  concept,  and  a  connection  between  the  semi-  and  fully-discrete  stability 
must  be  used.  The  third  theorem  provides  this  connection. 

Theorem  3:  Under  mild  restrictions  [9],  if  a  semi-discrete  approximation  is  stable  in  a 
generalized  sense  and  a  Runge-Kutta  method  which  is  locally  stable  is  used  to  time  march 
the  semi-discretization,  the  resulting  totally  discrete  approximation  is  stable  in  the  same 
sense  so  long  as  the  stability  region  of  the  R-K  method  encompasses  the  norm  of  the  semi¬ 
discretization. 

It  should  be  noted  [9]  that  the  stability  definitions  used  in  the  first  two  theorems  (G-K-S 
stability)  are  different  than  that  of  the  third  (generalized  stability).  The  first  two  theorems 
rely  on  G-K-S  stability  (sometimes  referred  to  as  the  Kreiss  condition)  in  the  semi-  or  fully- 
discrete  case.  At  least  two  different  definitions  of  G-K-S  stability  are  encountered  in  the 
cumulative  works  on  G-K-S  analysis.  They  are: 

Definition  1:  The  IBVP  is  stable  if  for  t]  >  ?/0,  the  solutions  of  eqn  (1)  with  homogeneous 
initial  data  satisfy 

l#(0,  S)f  +  (v  -  *?0)2||U(.,S)||2  <  I<(\g\2  +  \\F(.,S)\\2)  (16) 

where  t/0  and  K  are  universal  constants,  and  the* denotes  the  Laplace  transformed  variables. 

Definition  2:  The  IBVP  is  stable  if  for  t]  >  r/c,  the  solutions  of  eqn  (1)  with  homogeneous 
initial  data  and  g  =  0  satisfy 

(v-Vo)\\U(.,S)\\<K(\\F(.,S)\\).  (17) 

Both  of  these  definitions  can  be  related  to  what  is  often  referred  to  as  Lax  stability,  which 
is  the  numerical  equivalent  to  an  energy  estimate  satisfied  by  the  discrete  form  of  eqn  (1,2). 

Theorem  3  relies  on  a  different  definition  of  stability  that  can  be  expressed  as 

Definition  3:  The  IBVP  is  stable  if  for  homogeneous  boundary  data  (g=0)  an  estimate 
of  the  form 

\\U(;t)\\  =  Vfcxp«<‘— >(||C/(.,l.)||  +  /'  ||F(.,T)||dr).  (IS) 

Jo 

Definition  2  is  the  most  restrictive  of  the  three  conditions,  but  it  has  been  conjectured 
by  Kreiss  [9]  that  Definition  2  and  3  are  equivalent.  The  subtle  differences  between  these 
definitions  of  stability  will  not  affect  the  conclusions  in  this  work,  and  the  terms  G-K-S 
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stability  and  Lax-stability  shall  be  used  interchangeably.  We  now  describe  in  some  detail 
the  implications  of  the  three  Theorems. 

lhe  chief  difficulty  in  stability  analysis  is  not  which  definition  is  chosen,  but  rather 
applying  that  definition  to  practical  numerical  methods  and  obtaining  a  stability  bound. 
Since  Fourier  methods  are  not  well  suited  on  the  finite  domain,  stability  analysis  can  be 
carried  out  by  energy  methods  or  by  normal  mode  analysis.  Energy  methods  are,  in  general, 
very  difficult  to  perform  on  high-order  schemes.  The  modal  relationships  are  simple  to  define 
but  analytic  solutions  are  often  intractable. 

Theorem  1  describes  how  G-K-S  analysis  can  be  used  to  augment  finite  domain  modal 
analysis.  The  original  finite  domain  modal  analysis  is  broken  into  the  analysis  of  three 
equivalent,  yet  simpler  modal  problems.  Assuming  that  a  Cauchy  stable  scheme  is  used  for 
the  interior  grid-points,  the  inner-scheme  is  tested  for  stability  at  each  boundary  in  a  semi¬ 
infinite  spatial  domain.  In  so  doing,  the  stability  of  each  boundary  is  independent  of  the 
influence  from  the  other  boundary.  Stability  of  the  two  boundary  problems  implies  stability 
of  the  finite  domain  numerical  method.  In  addition,  Theorem  1  provides  a  “perturbation 
test”  for  “generalized  eigensolutions.”  The  test  establishes  the  stability  of  certain  borderline 
cases  in  rormal  mode  analysis. 

To  fully  appreciate  the  power  of  G-K-S  analysis,  a  normal  mode  analysis  of  the  fourth- 
order  compact  scheme  (4-4-4)  described  by  eqn  (13)  is  presented  for  the  coupled  finite  domain 
problem.  We  proceed  by  assuming  that  the  semi-discrete  problem  defined  in  eqn  (11)  has  a 
solution  of  the  form 


W)  =  exp51  h  (19) 

where  “S”  are  the  eigenvalues  of  the  matrix  M+  1  M% .  Substitution  into  eqn  (11)  yields  the 
generalized  eigenvalue  problem 

M+SVj  =  M+V  (20) 

(for  which  we  have  assumed  g(t)  =  0),  and  the  resolvent  equation  provided  by  the  inner- 
scheme  is 

(<t>j~\  +  4 (f>j  +  <f>j+i)S  =  3(— <f>j~\  +  <f>j+\)  j  =  2,  ...  ,iV  —  1  (21) 

where  S  =  ^-S.  It  is  apparent  that  <f>j  =  <j>0t will  satisfy  this  expression,  yielding  the 
following  equation  for  the  eigenvalues 

(-  +  4  +  k)S  =  3(-i  +  k).  (22) 

AC  AC 
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This  is  a  quadratic  expression  in  k,  and  there  will,  in  general,  be  two  solutions  which  will 
satisfy  eqn  (22).  The  eigenvectors  are,  therefore,  of  the  form  Oj  =  C\k^  +  C2k2j.  Simple 
manipulations  show  the  two  roots  are  related  by 

-2  -  K 

(23) 


K\  —  Kj  K'  2  — 


2  K  +  1 


and 


We  note  that  if 


4>j  =  C1K?  +  C2(-^rTT)- 


a  +  ib,  then 

|a'i  |  =  a 2  +  b‘ , 


*2 


2k  +  1 


b2  +  {a  +  2)2 
(2b)2  +  (2  a  +  2  )2 


(24) 


(25) 


and  | k j |  >  1  implies  |a2|  <  1  Thus,  each  a  is  dominant  near  one  of  the  boundaries  and  its 
influence  decays  monotonically  with  increasing  distance  from  that  boundary.  The  form  of  C\ 
and  C2  must  be  determined  from  the  boundary  conditions  at  j  =  l  and  j=N.  The  conditions 
at  the  two  boundaries  can  be  represented  as 


Ci[5-F1(a1)]  +  C'2[5-JF1(a2)]  =  0 
C,[5  -  FnM]  +  C2\S  -  Fn( a2)]  =  0 


(26) 


where  Fi  and  Fy  are  the  functional  relations  resulting  from  substituting  the  eigenvector 
and  eigenvalue  into  the  expressions  at  grid-points  j  =  l  and  j=N,  respectively.  Note  that 
S  =  S(k ),  and  that  every  term  is  a  function  of  k.  It  is  apparent  that  no  general  solution  to 
this  problem  exists  except  the  trivial  one.  Non-zero  solutions  exist  for  the  condition  where 
the  determinant  is  equal  to  zero.  The  determinant  condition  gives  the  resulting  expression 
for  k,  the  roots  of  which,  with  eqn  (22),  give  the  eigenvalues  and  eigenvectors  of  the  system. 

No  closed  form  expressions  are  known  for  the  roots  to  the  determinant  polynomial  in 
this  case  or  for  any  other  reasonable  boundary  closures.  The  numerical  solution  of  the 
determinant  polynomial  involves  finding  the  roots  of  an  “nlh”  order  polynomial  in  a.  It 
is  apparent  that  the  boundary  conditions  dramatically  influence  the  eigenvalue  spectrum 
in  the  matrix  M.  Only  in  the  limit  N  — >  oo  do  the  boundary  conditions  decouple.  The 
power  of  G-K-S  theory  comes  from  breaking  the  normal  mode  problem  into  three  separate 
problems.  The  roots  to  the  k  polynomials  do  not  depend  on  the  boundary  to  boundary 
coupling  prevalent  in  eqns  (25-26). 

Theorem  2  describes  what  constitutes  stability  for  the  two  IBVP’s  in  Theorem  1  for  the 
fully-  or  the  semi-discrete  case,  and  states  that  eqn(l)  must  satisfy  the  condition  that  no 
eigenvalues  or  generalized  eigenvalues  should  exist  for  7£(S)  >  r/0. 


Both  theorems  1  and  2  rely  on  a  definition  of  an  eigensolution  for  their  quarter-plane 
analysis.  Here,  an  eigensolution  is  presented  for  the  semi-discrete  case.  Similar  definitions 
exist  for  the  fully  discrete  case.  Referring  to  the  semi-discretization  defined  by  eqn  (5), 

Definition  f.  An  eigensolution  for  the  IBVP  defined  by  eqn  (5)  is  the  nontrivial  function 
V(x,s)  satisfying  [8]  : 

I  SV  =  M  V  x  >  0 
II  V//(0,s)  -  Tl/7/(0,s)  =  0 

III  K{S)  >  0 

IV  for  7 Z(S)  >0  ,  V(x,s)  is  bounded  as  x  — »  oo 

V  for  TZ  =0  and  |«|  =  1  ,  a  perturbation  inside  the  unit  circle  of  n  (|/c|  =  1  —  t,e  >0) 
cannot  produce  an  eigenvalue  71  >  8,8  >0. 

Note  that  since  ip  —  xp0KJ >  condition  (IV)  implies  |/c|  <  1. 

We  refer  to  an  eigensolution  of  the  form  (IV)  or  (V)  as  a  G-K-S  eigenvalue  or  a  generalized 
G-K-S  eigenvalue,  respectively.  With  these  conditions,  the  test  for  numerical  stability  has 
been  simplified  from  the  coupled  normal  mode  analysis  to  tests  involving  7Z(S)  >  0,  for 
|k|  <  1  at  each  boundary,  plus  the  exceptional  case  when  both  77.(5)  =  0,  and  |k|  = 
1.  Theorem  3  relates  the  stability  of  the  semi-discretization  to  the  stability  of  the  fully 
discrete  numerical  method.  Obviously,  analogous  stability  definitions  must  be  chosen  for 
both  the  semi-  and  fully-discrete  cases,  specifically  the  generalized  stability  condition  given 
by  Definition  3.  Theorem  3  relies  on  temporal  advancement  schemes  which  are  locally  stable 
numerical  methods.  For  a  locally  stable  numerical  method,  the  stability  domain  |z|  <  1 
(z  is  the  amplification  factor)  in  the  complex  plane,  encompasses  within  the  LH-P  an  open 
semi-circle  of  radius  “Rf'  centered  at  the  origin,  and  symmetric  about  the  real  axis.  The 
standard  fourth-order  R-K  method  satisfies  this  condition.  Discretizing  time  with  the  fourth- 
order  R-K  method  in  eqn  (5)  produces  a  fully  discrete  method  defined  by 

VJ(t  +  k)  =  L(kM)VJ(t)  +  L(k)BJg(t)]  j  =  1,  ...  ,JV  (27) 

where  the  time-step  At  =  k,  and  where  “L(kM)”  is  the  polynomial  in  “k  M”  describing 
the  time  discretization.  Then  under  very  mild  restrictions  (see  Kreiss  [9])  on  the  eigenvalue 
structure  of  the  matrix  “L(kM)”,  if  the  semi-discrete  approximation  is  stable  in  a  generalized 
sense,  the  totally  discrete  approximation  is  stable  in  the  same  sense  for  ||fcM||  <  R\. 
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We  have  outlined  a  systematic  approach  to  addressing  the  finite  domain  stability  problem 
for  the  fully-discrete  numerical  approximation  to  eqn  (1).  The  remainder  of  this  work  will 
describe  the  application  of  these  techniques  to  several  higher-order  finite  difference  schemes. 


5.  FOURTH-ORDER  BOUNDARY  CONDITIONS 


Before  proceeding  with  the  stability  analysis  of  various  higher-order  boundary  conditions, 
an  example  is  presented  that  illustrates  the  necessity  of  (N  —  l)(/i  order  boundary  closure 
for  an  Nth  order  inner-scheme.  The  example  also  provides  a  numerical  test  to  verify  that 
the  G-K-S  theory  is  accurately  predicting  the  stability  behavior  of  the  various  numerical 
schemes.  Consider  the  method-of  lines  approximation  to  the  scalar  wave  equation 


dU_  dU_ 
Dt  dx 


=  0,  -1  <  X  <  1,  t  >  0 


(28) 


U(t,  —  1)  =  sin  2tt(  —  1  —  t),  U(0,  x)  =  sin  2ttx  —  1  <  x  <  1,  f>0  (29) 

where  the  spatial  discretization  is  accomplished  by  the  fourth-order  compact  scheme  de¬ 
scribed  in  detail  by  eqn  (12).  The  exact  solution  is 


U{t,x)  =  sin27r(x  —  f),  —  1  <  x  <  1,  t  >  0. 


(30) 


The  comparison  of  the  exact  solutions  with  that  obtained  from  various  lower-order  closures 
to  the  fourth-order  inner-scheme  provides  a  measure  of  the  boundary  influence.  We  begin 
by  showing  a  grid  convergence  study  performed  to  show  the  formal  accuracy  of  the  resulting 
schemes.  The  boundary  condition  formulas  expressed  at  the  grid-point  “j=0”  were 


d Vo  ^ 0 V\  —  5 Vo  -f-  4 Vj  +  V‘2 
dx  ~  dx  2Ax 

dV0  d\\  -V0  +  V7! 

dx  dx  2Ax 

d  Vo  —  Vo  +  Vi 

dx  Ax 


(31) 

(32) 

(33) 


which  represent  third-,  second-,  and  first-order  closure  at  the  inflow  boundary,  respectively. 
As  mentioned  earlier,  the  physical  boundary  condition  was  imposed  at  the  point  tkj=0,”  and 
the  actual  closure  occurs  at  the  point  “j  =  T’.  The  closure  could  have  been  written  explicitly 
for  the  point  “j  =  l”  by  combining  these  formulas  with  the  inner-scheme  at  “j  =  l” 


dy0  d\\  c/W  _  -3  Vp  +  3  Vj 
Ox  ^  dx  +  dx  Ax 


(34) 


13 


resulting  in  formula  of  the  form 


d\\ 

dv2 

-V  'o  -  4  V  A  +  5  V4 

9 

~  dx 

+ 

dx 

2  Ax 

(35) 

dVi 

dv2 

—  V'o  —  2  V  i  T  3V2 

(36) 

dx 

+ 

dx 

Ax 

,  d\\ 

dV2 

—2  V0  -  Vt  +  31  g 

(37) 

4  dx 

+ 

dx 

Ax 

At  the  outflow  boundary,  closure  was  accomplished  with  the  expressions 


dd/V  ^dlbv-i 
dx  “  dx 

ovN  dvN-i 

3x  dx 

dvN 

dx 


— 5V,v  +  4VV-1  T  Vy -  2 
%Ax 

—  l/v  +  V/v-i 

2Ax 

—  V,v  +  V/v-i 

Ar 


(38) 

(39) 

(40) 


which  represent  third-,  second-,  and  first-order  spatial  accuracy,  respectively.  These  expres¬ 
sions  are  valid  for  the  point  "j=N”  since  there  is  no  physical  condition  to  impose  there. 

In  all  cases,  the  temporal  discretization  was  accomplished  with  the  fourth-order  Runge- 
Kutta  algorithm.  At  every  iteration,  the  solution  was  advanced  for  the  grid-points  “j=0,N”, 
followed  by  the  imposition  of  the  boundary  condition  at  the  point  “j=0.”  Since  the  inflow 
boundary  condition  was  a  nonlinear  function  of  time,  care  was  taken  to  evaluate  the  function 
corresponding  to  the  proper  intermediate  level  in  time.  Failure  to  do  so  degraded  the  formal 
accuracy  of  the  method.  The  CFL  used  in  the  simulations  was  in  the  range  0.1  <  CFL  <  1. 
In  no  case  did  this  violate  the  Von  Neumann  stability  condition  for  the  Cauchy  problem.  It 
should  be  noted  that  the  formal  truncation  of  the  method  is  0(A(4),  and  the  error  in  time 
decays  to  the  fourth  power  of  the  CFL  for  a  given  grid,  and  can  be  made  as  small  as  desired 
by  decreasing  the  CFL.  Typically,  CFL’s  <0.1  were  not  needed  since  the  dominant  terms 
in  the  modified  equation  were  negligible  compared  with  spatial  terms.  Further  reduction  of 
the  CFL  resulted  in  no  change  in  the  error  of  the  scheme.  Recognition  of  this  fact  allows  one 
to  test  the  formal  accuracy  of  methods  with  spatial  accuracy  higher  than  fourth-order  for 
sufficiently  small  values  of  the  CFL,  and  was  used  in  the  sixth-order  simulations  to  determine 
formal  accuracy.  Finally,  a  third-order  Runge-Kutta  was  used  to  test  the  generality  of  the 
temporal  discretization  in  several  cases.  The  results  were  quantitatively  similar. 

The  simulations  were  all  run  to  equivalent  times  T— 25,  for  all  grids  and  methods  at  a 
CFL  of  0.25.  The  error  at  T  was  then  calculated  and  reported  as  an  L2  norm.  The  L <*>  norm 
produced  similar  results,  but  is  not  reported  here.  For  methods  which  are  Lax-stable,  the 
error  is  bounded  uniformly  at  each  stage  for  0  <  A/  <  r  and  0  <  kAt  <  T,  where  k  is  the 
number  of  time-steps.  Doubling  the  grid  at  constant  CFL,  should  decrease  the  error  at  time 
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level  T  by  a  factor  |p  where  "p '  is  the  order  ol  the  method.  The  formal  accuracy  ot  each 
scheme  was  determined  tor  each  of  the  Lax-stable  schemes  in  this  manner. 

Figures  1-3  show  the  results  of  this  study.  The  log10  of  the  L2  error  is  plotted  as  a 
function  ot  the  log10  ot  the  number  of  grid-points  resolving  one  period  of  the  sine  wave. 
The  grid  density  ranges  from  10  to  25  grid-points/(2/r  radians).  Note  that  once  a  threshold 
accuracy  is  attained  th  oints  appear  to  decrease  linearly.  The  slope  of  the  data  (dY/dX) 
gives  the  apparent  formal  accuracy  ot  the  method.  In  all  cases,  formal  fourth-order  accuracy 
is  obtained  with  third-  or  tourth-order  boundary  conditions.  Second  order  closure  results 
in  third-order  formal  accuracy,  and  first-order  closure  results  in  second-order  accuracy.  It 
is  also  apparent  that  the  accuracy  is  greater  for  the  fourth-order  than  for  the  third-order 
closure,  regardless  ot  the  formal  accuracy.  These  results  are  in  agreement  with  the  theory  of 
Gustafsson  [3]  predicting  that  to  retain  formal  accuracy  of  a  numerical  method,  boundary 
conditions  of  order  (N-l)  must  be  imposed  to  retain  \nh  order.  Another  interesting  feature  is 
that  imposition  of  lower-order  boundary  conditions  at  the  outflow  plane  results  in  a  greater 
degree  of  error  than  at  the  inflow  plane.  One  might  have  thought  that  since  the  characteristic 
is  pointing  out  of  the  domain  at  the  outflow  boundary,  error  would  be  swept  immediately 
out  of  the  domain.  This  is  obviously  not  the  case. 

We  will  now  derive  the  formal  stability  of  the  numerical  boundary  conditions  used  in 
this  example.  It  can  be  shown  that  the  fourth-order  compact  scheme  is  “Cauchy  stable” 
for  C’FL  <  1.63.  The  stability  of  the  inflow  and  outflow  boundary  conditions  on  the  semi¬ 
infinite  domain  must  be  demonstrated.  We  begin  by  testing  the  outflow  stability.  The  partial 
differential  equation  is 

DU  dU 

-—-  —  =  Q,  x  >  0,  t>  0.  (41) 

dt  dx  ’ 

No  boundary  condition  is  required  in  this  problem,  although  a  NBS  is  imposed  at  x  =  0.  As 
was  done  on  the  normal  mode  analysis  for  the  finite  domain,  we  assume  a  solution  of  the 
form  Vj{t)  =  expM0j,  where  d)j  =  <f>0KJ .  Substitution  into  the  inner  scheme  produces  the 
resolvent  condition  for  the  eigenvalue  S 

(-  +  4  +  k)5  =  3(—  +  k)  (42) 

K  K 

at  each  grid-point  j  >  1.  At  grid-point  j=0,  the  scheme  was  closed  with  one  of  the  bound¬ 
ary  expressions  given  in  eqns  (12,  38,  39,  40)  (obviously  written  for  the  grid-point  j=0). 
Substitution  of  Vj(f)  into  these  expressions  yields 


-17  +  9k  +  9s2  -  x-3 

-  ■ )  +  4  k  +  K 


(6+  18k). S’ 

(2  +  4k)N 


■) 


(43) 

(44) 


-2  +  2k 

—  1  +  K . 


(1  +  K  )S  =  -2  +  2  K  (45) 

5  =  —1  +  K.  (46) 

Solving  eqn  (42)  for  5  and  substituting  into  eqns  (43,  44,  45,  46)  yield  polynomials  in  k  of 
the  form 

(4+4+1)  =°  (4“)  (47» 

=  °  <3">  <48> 

(k-1)3 

(,'+4.4l)=0  <2  )  ‘49> 

(k  —  1)2(k  +  2) 

W4+r=° (1  »■  (50) 

Clearly,  the  only  value  of  k  which  will  simultaneously  satisfy  the  inner  resolvent  condition 

and  any  of  the  outflow  boundary  conditions  is  k  =  1.  Substitutions  of  k  =  1  into  the 

resolvent  condition  produces  5  =  0,  the  case  for  which  the  perturbation  test  in  G-K-S 
theory  (condition  V)  must  be  used  to  show  stability.  Substituting  k  —  1  —  t  into  the 
resolvent  condition  produces  to  first  order 

65  =  — 6e.  (51) 

For  e  >  0  and  |k|  <  1,  5  <  0  showing  stability  of  the  perturbation.  All  the  tested  outflow 
boundary  conditions  are  G-K-S  stable.  Thus,  k  =  1  is  not  a  generalized  eigenvalue. 

To  show  stability  of  the  inflow  conditions  we  study  the  partial  differential  equation 


=  0  (ls(). 


dU  dU  n 
dt  +  dx  ~  °’  x  -  °’  1  -  0; 
with  the  boundary  condition  imposed  at  x  =  0  of  the  form 

V{0,t)=g(t).  t>  0.  (53) 

In  spite  of  the  physical  boundary  condition  being  imposed  at  j =0,  a  NBS  must  be  imposed 
at  j  =  1 ,  and  must  be  tested  for  stability.  Substituting  Vj(t)  =  exp5t  <j>j,  where  <f>j  —  <f>0K* ,  into 
the  inner-scheme  produces  the  resolvent  condition  identical  to  eqn  (42).  Substitution  into 
the  first  through  fourth-order  boundary  conditions  defined  by  eqns  (12,  35,  36,  37),  produce 
equations  for  5  of  the  form 


(6  +  6k)5  = 

-9  +  9k  +  k2 

(54) 

(4  +  2k)5  = 

—4  +  5k 

(55) 

(6  +  k)5  = 

—4  +  6k 

(56) 

(4  +  2k)  5  = 

—2  +  6k. 

(57) 

Note  that  without  loss  of  generality,  we  have  assumed  g(t)  =  0  on  the  inflow  boundary  j=0. 
This  eliminates  the  influence  of  j=0  in  the  boundary  polynomial,  and  reduces  the  order  of 
the  polynomials  by  one.  Solving  the  resolvent  equation  for  S  and  substituting  into  eqns  (54, 
55,  56,  57)  yields  polynomials  in  k  of  the  form 

(a-4  -  5k3  +  10k2  -  9k  +  9) 

(k2  +  4k  -f  1) 

(k3  —  4k2  +  5k  —  8) 

(k2  -j-  4k  +  1 ) 

(k2  —  2k  +  7) 

(k2  +  4k  +  1 ) 

(k2  —  2k  —  1 1) 

(k2  +  4k  +  1 ) 


the  roots  of  which  are 


k  =  2.286  ±  1.215*,  0.3134  ±  1.1 38i 

(4"*) 

(63) 

k  =  3.218,0.3906  ±  1.527* 

(3rJ) 

(64) 

k  =  1  ±  v/6  i 

(2nd) 

(65) 

K  =  1  ±  2 \/Si 

(1st) 

(66) 

=  0 

(4'*) 

(53) 

=  0 

(3'J) 

(59) 

=  0 

(2nd) 

(60) 

=  0 

(1") 

(61) 

where  i  =  \/—  l .  It  is  apparent  that  |k|  >  1  in  all  of  these  expressions.  There  are  no 
eigenvalues  or  generalized  eigenvalues,  and  the  inflow  boundary  is  stable  to  these  closures. 
Given  that  the  Cauchy  problem  and  the  two  quarter-plane  problems  are  stable,  implies  that 
the  finite  domain  problem  defined  in  eqn  (12)  is  G-K-S  stable  for  all  boundary  conditions 
specified  thus  far. 

A  more  rigorous  example  of  the  ability  of  the  G-K-S  theory  to  predict  the  stability  of 
the  fourth-order  compact  scheme  is  demonstrated  by  a  pathological  inflow  boundary  scheme. 
Since  the  inflow  problem  involves  a  NBS  at  the  grid-point  j  =  l  which  is  biased  in  the  “down¬ 
wind”  direction,  we  suspect  that  it  will  be  more  sensitive  to  instability  than  the  outflow 
boundary.  We  formulate  a  boundary  scheme,  which  is  a  linear  combination  of  the  first- 
and  second-order  schemes,  noting  that  they  both  have  been  shown  to  be  stable  boundary 
treatments.  The  resulting  scheme  which  we  shall  use  at  the  inflow  boundary  is 


(1+2„^  +  ^  =  2,1 


dx  dx 


Ax 


(67) 


where  f3  is  the  amount  of  first-order  influence  in  the  formula.  For  (3  —  0,  the  standard  second- 
order  closure  is  obtained.  For  all  other  values  of  /?,  the  scheme  is  formally  first-order  accurate. 
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We  shall  denote  this  scheme  as  ( 1 1  -4-4)  since  the  inflow  boundary  is  a  one-parameter  family 
of  first  order  schemes,  and  the  outflow  is  closed  with  a  fourth-order  formula.  As  was  done 
previously,  the  dependence  of  the  formula  on  the  space  derivative  at  j=0  is  eliminated  from 
eqn  (67)  by  the  inner-scheme  written  at  j  =  l.  The  resulting  expression  is 

ri,,,  ,  hdVt  (1+2 l3)dV2  — (1  +  4/? ) Vo  —  2(1  +  3)V\  +  3(1  +  2/?) /cox 

l2(1  +  2J)  2|37T— 2~ &T  =  - 2^ - '(68) 

Substituting  V}(t]  —  expM</q,  </q  =  4>0k]  into  eqn  (68)  and  noting  the  expression  must  be 
valid  for  the  case  g(t)  =  0,  results  in  a  boundary  scheme  of  the  form 

([2(1  +  2 0)  -  i]  +  — =  -(1  +  D)  +  |(1  +  2/3)*.  (69) 

Solving  for  5  from  the  inner  scheme  eqn  (42)  and  substituting  it  into  eqn  (69)  produces  a 
polynomial  in  k  of  the  form 

(-1  +  2 i3)k2  +  (2  -  4 0)k  -  (7  +  22/?)  =  0.  (70) 


Solving  for  k  yields 


k  =  l  ±y/6 


\ 


(4/?  +1)  1 

(2jTTy 


(71) 


A  double  root  exists  for  f3  =  ^  and  k  =  1.  Substituting  k  =  1  into  the  resolvent  expres¬ 
sion  yields  S  =  0,  and  a  perturbation  test  shows  that  the  boundary  exhibits  a  generalized 
eigenvalue  instability.  Further  inspection  of  eqn  (71)  shows  that  |«|  <  1  over  the  range 
“  <  /?  <  All  other  values  of  j3  need  not  be  considered  as  candidates  for  instability  since 
|k|  >  1.  Substitution  of  the  expression  for  k  obtained  from  eqn  (71),  into  eqn  (69),  yields 
an  expression  for  S  in  terms  of  the  parameter  j3.  Numerical  evaluation  of  this  expression 
shows  that  1Z(S)  >  0  for  —.37  <  /?  <  Thus,  an  eigensolution  exists  for  this  range,  and 
the  coupled  inner-boundary  scheme  is  unstable. 

To  verify  these  findings,  the  model  scalar  wave  equation  described  by  eqns  (28,  29,  30), 
was  solved  with  the  pathological  inflow  boundary  conditions  describe  in  eqn  (67).  Fourth- 
order  boundary  conditions  were  used  at  the  outflow  boundary.  Figure  (4)  shows  the  results  of 
the  numerical  investigation.  Plotted  is  the  log10  of  the  L2  error  of  the  solution  integrated  to 
a  fixed  time  "T,”  as  a  function  of  the  parameter  /?,  ranging  from  [+p,  \).  Two  grid  densities 
are  shown  in  the  study,  and  behave  similarly.  The  theoretically  predicted  range  of  instability 
—  .37  <  /?  <  —  is  replicated  in  the  numerical  study  to  within  graphical  limitations  of  the 
plot.  In  the  unstable  regime,  the  error  grew  exponentially  with  the  number  of  iterations 
required  to  reach  the  time  level  T,  and  quickly  became  very  large. 
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Although  those  boundary  conditions  are  pathological,  this  study  points  to  the  fact  that 
imposition  of  a  first-order  inflow  boundary  conditions  is  not  sufficient  to  guarantee  stability 
with  the  fourth-order  compact  inner-scheme.  A  similar  experiment  was  performed  on  the 
outflow  boundary  using  a  linear  combination  of  first-  and  second-order  boundary  conditions. 
For  those  boundary  conditions,  no  eigensolutions  could  be  found.  That  is  not  to  imply  that 
a  completely  arbitrary  outflow  boundary  condition  of  at  least  first-order  accuracy  is  always 
stable,  ft  does  indicate  that  the  outflow  is  less  susceptible  to  instabilities  than  the  inflow 
boundary. 

Another  note  is  appropriate,  concerning  the  complexity  of  the  G-K-S  analysis  for  compact 
schemes.  For  the  fourth-order  compact  inner-scheme,  the  polynomial  of  the  highest  degree 
which  could  not  be  factored  (resulting  in  a  numerical  solution)  was  of  degree  four,  and 
resulted  from  the  closure  with  the  fourth-order  boundary  conditions  at  the  inflow  plane.  The 
G-K-S  analysis  of  the  fourth-order  explicit  scheme  described  by  eqn  (9)  was  performed  by 
Stiikv.erdc*  [3].  Stability  polynomials  of  order  eight  were  obtained  for  fourth-order  closure 
at  the  inflow  plane.  It  is  apparent  that  the  stability  polynomials  resulting  from  compact 
stencils  are  simpler  expressions.  This  simplicity  (and  MACSYMA)  will  allow  us  to  analyze 
schemes  of  sixth-order  spatial  accuracy  without  insurmountable  algebraic  polynomials  at 
each  boundary. 

6.  SEMI-DISCRETE  EIGENVALUE  ANALYSIS 

The  uniformly  fourth-order  explicit  scheme  (4, 4-4-4, 4)  analyzed  by  Strikwerda  [8]  and  the 
uniformly  fourth-order  compact  scheme  (4-4-4)  presented  here  are  both  G-K-S  stable  for  the 
semi-discr^e  problem  and,  therefore,  will  exhibit  generalized  stability  for  the  fully  discrete 
problem  if  advanced  with  a  locally  stable  temporal  scheme.  This  definition  of  stability  ensures 
that  the  error  of  the  numerical  solution  will  remain  uniformly  bounded  for  all  times  by  an 
exponentially  growing  amount.  The  exponential  growth  rate  of  the  error  is  asymptotically  (N 
— >  oo,  where  N  is  the  total  number  of  grid-points  used)  independent  of  the  grid  used.  Thus, 
grid  refinement  studies  with  these  methods,  performed  by  integrating  the  governing  equation 
to  a  fixed  time  level  T  on  successively  finer  grids  will  demonstrate  that  the  numerical  solution 
converges  to  the  exact  solution  at  a  rate  of  at  least  the  order  of  the  m’ethod. 

A  disturbing  feature  of  this  stability  definition  is  that  the  solution  is  not  required  to 
remain  bounded  for  all  times,  even  though  the  physical  solution  remains  bounded  for  all 
times.  Figure  (5)  shows  a  grid  refinement  study  performed  with  the  fourth-order  compact 
(4-4-4)  scheme  demonstrating  this  behavior.  The  model  equation  was  that  described  by 
eqns  (28,  29,  30);  the  time  interval  was  0  <  t  <  100  and  the  grids  used  were  21,  41,  and  81 
grid-points,  respectively.  Time  was  advanced  with  a  fourth-order  Runge-Kutta  scheme  in  all 
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cases.  The  exact  solution  is  a  traveling  sine  wave  of  amplitude  1,  for  all  times.  Shown  is  the 
log  10  of  the  L 2  error,  plotted  as  a  function  of  time.  Simulations  on  all  three  grids  were  run 
at  various  CFL’s.  The  initial  portion  of  the  simulation  is  characterized  by  nearly  constant 
levels  of  error  on  all  three  grids.  After  a  sufficiently  long  time,  the  unstable  modes  in  the 
numerical  solution  dominate  the  spatial  truncation  error.  From  that  point  on,  the  solution 
diverges  exponentially  from  the  exact  solution.  Note  that  the  growth  rate  in  time,  of  the 
unstable  modes  of  the  solution  is  nearly  grid  independent,  and  that  at  any  time  T  refining 
the  grid  by  a  factor  of  two  results  in  a  factor  of  16  decrease  in  the  error.  It  is  also  evident 
that  at  large  times  the  actual  error  will  be  exponentially  large.  An  interesting  feature  of 
the  numerical  method  is  that  the  exponential  growth  of  the  solution  is  dependent  on  the 
CFL  used  to  advance  the  solution.  For  CFL  =  1,  the  solution  does  not  grow  in  time,  while 
for  CFL  <  a  (a  is  some  critical  value  less  than  CFLmar),  exponential  growth  is  observed. 
(This  feature  will  be  explained  later,  in  terms  of  the  amplification  factor  of  the  scheme.) 
Regardless  of  the  CFL,  fourth-order  convergence  is  observed  with  the  scheme. 

To  understand  the  fundamental  nature  of  the  fixed  grid,  temporal  divergence  of  the 
solution  in  the  previous  example,  it  is  instructive  to  study  the  eigenvalue  spectrum  of  the 
spatial  discretization  operator.  As  a  semi-discretization,  eqn  (28)  can  be  written  in  the  form 
of  eqn  (5)  as 

f  =  Ml' +  J  =  l,  ...,N  (72) 

where  M  is  the  N  x  N  matrix  describing  the  spatial  discretization  operator,  and  Bj  g(t) 
represents  the  physical  boundary  data.  Assuming  that 

P~1MP  =  S ;  P~lU  =  V\  P~1Bg(t)  =  H  (73) 

where  S  is  a  diagonal  matrix,  and  P-1  and  P  are  similarity  transforms  composed  of  the  left 
and  right  eigenvectors  of  the  matrix  M,  respectively,  eqn  (72)  takes  the  form 

—  SUj  +  H:,  j  =  1,  ...,iV.  (74) 

The  solution  to  eqn  (74)  is 

Uj{t)  =  expsT  Uj{ 0)  +  /  exps;(t_r)  //j(t)</t,  j  =  1,  ... ,  N.  (75) 

Jo 

In  this  form,  the  solution  to  eqn  (74)  depends  exponentially  on  the  eigenvalues  Sj  of  the 
matrix  M.  Note  that  this  solution  assumes  that  the  eigenvalues  S}  are  not  degenerate,  and 
that  H(t)  is  not  at  a  resonance  frequency.  If  either  of  these  situations  occur,  then  the  solution 
would  include  terms  proportional  to  tp  ex pS}t,  where  “p”  is  the  order  of  the  degeneracy.  The 
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precise  behavior  depends  on  the  temporal  nature  ot  II ( t ) ,  but  in  general,  for  boundary  data 
which  remains  bounder!  for  all  time,  the  solution  grows  only  for  the  modes  which  have 
eigenvalues  ,Sj  with  positive  real  parts.  In  addition,  the  growth  rate  will  be  governed  by 
the  eigenvalue  with  the  largest  positive  component.  Thus,  any  spatial  discretization  to  the 
semi-discrete  problem  defined  in  eqn  (72)  will  exhibit  exponential  divergence  of  the  solution 
from  the  bounded  physical  solution,  if  it  has  an  eigenvalue  in  the  right  half  of  the  complex 
plane. 

Figures  (6-7)  show  the  eigenvalue  spectrum  resulting  from  the  explicit  fourth-order  and 
the  compact  fourth-order  spatial  operators,  closed  at  the  boundaries  with  schemes  of  third- 
or  fourth-order  accuracy.  In  shorthand  nomenclature,  the  explicit  cases  (3,3- 4-3,3)  and  (4,4- 
4-4,4)  are  shown  in  Figure  (6),  and  the  compact  cases  (3-4-3)  and  (4-4-4)  are  shown  in 
Figure  (7).  The  spectrums  are  shown  on  grids  of  21,  41,  and  Si  points,  respectively.  Note 
that  closing  the  inner-schemes  with  third  order  MBS  in  both  cases  results  in  an  eigenvalue 
spectrum  which  is  bounded  to  the  Left  Half-Plane  (LH-P),  and  that  the  uniformly  fourth- 
order  schemes  cross  over  the  imaginary  axis  into  the  Right  Half-Plane  (RH-P)  of  the  complex 
plane. 

For  long  times,  the  maximum  eigenvalues  from  the  uniformly  fourth-order  schemes  very 
accurately  predict  the  exponential  growth  of  the  solution.  In  Figure  (5),  the  solutions  ob¬ 
tained  from  the  (4-4-4)  compact  scheme  grow  exponentially  in  time.  Assuming  the  er¬ 
ror  can  be  represented  functionally  as  c,v(0  =  t  ,v(0)  expQ/v(,  where  N  is  the  number  of 
grid-points  used  in  the  spatial  discretization,  a  growth  rate  a,\  can  be  determined  numer¬ 
ically  Similarly,  from  an  eigenvalue  determination,  an  effective  grow  rate  as  defined  by 
exp“'‘VAf  =  |6'm,1J.(A()|A ,  can  be  calculated,  where  Gmax  is  the  numerical  amplification  ob¬ 
tained  from  the  temporal  advancement  scheme.  For  the  fourth-order  Runge-Kutta  scheme 


Ci  j  —  1  +  Sj  CSt  + 


(SjAt2)  (SjAts)  (Sj&t 


9! 


+ 


3! 


+ 


J  =  1, 


(76) 


and  \G max |  will  frequently  correspond  to  the  maximum  eigenvalue  'R{Smax)-  Table  (1)  shows 
a  comparison  of  the  observed  growth  rate  of  the  (4-4-4)  compact  scheme  with  that  predicted 
from  an  eigenvalue  determination.  In  each  case,  the  maximum  eigenvalue  is  used  to  predict 
the  temporal  grow  of  the  solution. 


Grid 

ume.rical 

°K(.S„,ar) 

21 

0.1321 

0.1315 

41 

0.1476 

0.1474 

81 

0.1537 

0.1479 

Table  1:  Numerical  vs.  Theoretical  Growth  Rate;  (4-4-4). 
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The  agreement  is  very  good,  with  a  slight  discrepancy  in  the  comparison  on  the  81  grid-point 
case.  In  Figure  (5),  note  the  oscillatory  nature  of  the  growth  of  the  solution.  The  uncertainty 
of  the  phase  of  the  solution  accounts  for  the  discrepancy  in  the  predicted  growth  rate  in  that 
case. 

A  necessary  condition  for  Lax  stability  of  the  finite  domain  semi-discretization  can  be 
expressed  in  terms  of  the  eigenvalues  of  the  spatial  matrix  operator  as 

^>0;  i  =  l,  ...,JV  (77) 

where  Sj  are  the  eigenvalues  of  the  spatial  operator,  and  N  is  an  arbitrary  number.  It  is 
apparent  that  the  eigenvalue  structure  asymptotically  approaches  a  bound  in  the  RH-P  as 
N  — *  oo.  All  the  fourth-order  schemes  presented  thus  far  have  satisfied  this  constraint.  For 
the  third-order  NBS’s  (3, 3-4-3, 3)  and  (3-4-3),  the  constant  is  a  —  0,  while  for  the  fourth-order 
NBS’s,  the  constant  is  greater  than  zero. 

As  mentioned  earlier,  a  curious  feature  of  the  (4-4-4)  (as  well  as  other  high-order  spatial 
schemes)  is  that  the  growth  of  the  solution  is  CFL  dependent  on  all  grids.  For  CFL’s  close  to 
the  CFLmax ,  as  determined  from  Von  Neumann  stability  analysis,  the  schemes  are  bounded 
in  time.  For  sufficiently  small  CFL’s,  the  schemes  begin  to  diverge  exponentially  in  time. 
(Again,  it  should  be  emphasized  that  in  either  case,  the  scheme  is  still  G-K-S  or  Lax  stable.) 
This  behavior  can  be  explained  by  noting  a  particular  feature  of  the  fourth-order  Runge- 
Kutta  time  advancement  scheme  as  well  as  some  of  the  other  locally  stable  time  schemes. 
The  stability  bound  of  a  time  advancement  scheme  is  defined  as  the  locus  of  points  in  the 
complex  plane  where  \z\  <  1.  Clearly,  \z\  =  1  divides  the  plane  into  two  regions.  When 
the  spatial  eigenvalues  (scaled  by  At)  of  a  particular  discretization  lie  entirely  within  the 
\Z\  =  1  boundary,  the  combined  time-space  scheme  is  generally  stable.  The  stability  regime 
of  these  schemes  includes  a  semi-circular  portion  of  the  complex  plane,  centered  at  the  origin 
and  symmetric  about  the  real  axis  and  extending  into  the  left  half-plane.  In  addition,  they 
contain  a  small  part  of  the  right  half-plane,  although  not  near  the  origin.  If  the  spatial 
eigenvalues  which  lie  in  the  RH-P  are  encompassed  by  the  \Z\  =  1  boundary,  the  resulting 
scheme  is  stable.  If  a  At  is  chosen  such  that  the  \Z\  =  1  line  does  not  contain  the  RH-P 
eigenvalues,  then  the  solution  diverges  with  time.  Figure  (8)  shows  this  feature  for  the  (4-4- 
4)  spatial  scheme  and  the  fourth-order  Runge-Kutta  scheme  in  time.  For  CFL’s  which  are 
near  the  CFLmax,  the  maximum  amplification  rate  |Crmax|  is  less  than  one.  For  sufficiently 
small  CFL’s,  the  |(7max|  >  1  by  an  amount  that  is  proportional  to  7^.(5m0x),  and  the  solution 
will  diverge  exponentially  in  time.  It  can  be  shown  that  a  spatial  scheme  that  has  RH-P 
eigenvalues  can  always  be  made  to  diverge  exponentially  for  sufficiently  small  CFL’s  if  a 
conventional  third-  or  fourth-order  R-K  time  advancement  scheme  is  used. 
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7.  ASYMPTOTIC  STABILITY 


Lhe  previous  discussion  of  Figure  (5-8)  brings  out  a  subtle  point  in  the  Lax  stability 
theory.  In  dealing  vvitli  the  numerical  integration  of  time-dependent  partial  differential  equa¬ 
tions.  there  are  different  limit  processes  that  can  be  considered.  One  limit  is  the  behavior  of 
the  numerical  solution  as  the  mesh  size  Ax  — >  0  for  a  fixed  time  T.  Another  is  the  behavior 
of  the  solution  for  a  fixed  mesh  as  the  time  T  tends  to  infinity. 

Stability,  in  the  sense  of  Lax,  addresses  the  first  issue:  boundedness  of  the  numerical 
solutions  as  the  mesh  is  refined  at  a  fixed  physical  lime.  The  essence  of  the  Lax  equivalence 


theorem  is  that  if  the  numerical  solution  is  bounded  in  this  sense,  then  it  converges  to  the 
true  solution  in  the  limit  Ax  — »  0.  To  obtain  an  approximation  to  the  true  solution  at  time 
T,  one  integrates  the  IBVP  up  to  time  T  on  a  sequence  of  grids  as  Ax  — >  0.  This  sequence 
converges  to  the  exact  solution  for  all  time  levels  T. 

Nothing  in  this  definition  excludes  growth  in  time ,  and  specifically  allows  exponential 
growth  in  time  (see  eqn  (18)).  Moreover,  even  if  each  of  the  quarter-plane  problems  is  stable 
and  allow  no  growth,  in  time,  the  combined  finite  interval  problem  still  allows  exponential 
growth  in  time.  (The  Laplace  transforms  used  in  the  G-K-S  theory  are  legitimate  only  if 
growth  in  time  is  allowed.) 

Unfortunately,  for  genuinely  time  dependent  problems,  this  stability  definition  might  be 
too  weak,  in  particular  if  long  time  integration  is  being  carried  out.  The  reason  is  that  in 
order  to  achieve  any  reasonable  accuracy  for  large  times,  one  needs  an  excessive  number  of 
grid  points.  For  long  time  numerical  simulations  to  be  useful,  the  solution  of  the  semi-discrete 
problem  defined  in  eqn  (5),  must  be  bounded  in  time  as  well.  This  means  that  for  a  fixed 
mesh  N,  the  eigenvalues  of  the  matrix  M  in  eqn  (5)  have  non-positive  real  part,  and  those  with 
zero  real  part  have  (geometrical)  multiplicity  of  1.  This  is  called  Asymptotic  Stability. 

An  irksome  feature  of  asymptotic  stability  is  that,  by  itself,  asymptotic  stability  does  not 
imply  Lax  stability.  There  are  numerous  examples  in  the  literature  of  fully  discrete  schemes 
that  are  asymptotically  stable,  but  not  Lax  stable.  The  classic  example  of  this  is  the  case  of  a 
first  order  upwind  spatial  operator  being  advanced  with  an  Euler  explicit  time  advancement 
scheme.  The  eigenvalues  for  the  fully  discrete  system  are  (1  —  A),  occurring  a  degenerate  N 
times.  Eigenvalue  determination  suggests  that  the  CFL  of  the  scheme  should  be  2,  whereas. 
Von  Neumann  analysis  and  practical  experience  suggests  that  a  CFL  of  one  is  the  maximum 
stable  CFL.  The  discrepancy  can  be  explained  by  the  fact  that  for  1  <  CFL  <  2,  the  matrix 
norm  first  grows  rapidly  before  decaying  asymptotically  to  zero.  One  might  have  suspected 
this  difficulty  by  noting  that  in  the  semi-discrete  case,  the  degenerate  eigenvalues  give  rise  to 
geometric  growth  in  time,  only  later  to  be  dominated  by  the  exponentially  decaying  terms 
in  the  expressions. 
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An  example  for  the  semi-discrete  case  is  now  presented.  Consider  eqn  (5)  with  g(t)  =  0. 
Let  the  matrix  M  be  defined  by  the  second-order  central  difference  operator  for  the  interior 
points.  In  matrix  form,  M  can  be  written  as 


-1  1 

-1  0  1 


0 


-1  0  1 
-1  1 


(78) 


where  the  boundary  closures  at  grid-points  j  =  l  and  j=N  are  artificially  chosen.  Using  semi¬ 
discrete  modal  analysis,  we  determine  a  solution  to  eqn  (5)  by  assuming  the  form  Vj[t)  = 
exp St  4>j,  with  <f>  =  +  B{~)3 .  The  resolvent  condition  from  the  interior  scheme  yields 

Sj  =  |(/c  —  L).  Using  the  boundary  conditions  at  point  j=l  and  j=N  to  determine  the  values 
of  A  and  B  yields  the  expression  for  k  of  the  form 


(k  +  l)((-lf  -  k2")  =  0.  (79) 

The  roots  to  eqn  (79)  are  k  =  —  1  and  k  =  i  exp-^  for  j=l,N-l.  The  eigenvalues  are  thus: 
S  =  0  and  S  =  i  cos  jj-  for  j=l,N-l,  and  are  purely  imaginary.  The  spatial  discretization 
satisfies  our  definition  of  asymptotic  stability  for  values  of  N,  which  are  odd. 

The  spatial  discretization  defined  by  eqn  (78)  admits  a  generalized  eigenvalue  instability 
at  the  inflow  boundary.  Using  G-K-S  analysis  for  the  inflow  boundary  produces  compatibility 
equations  of  the  form 


25  k  =  k  —  1,  j  =  1 

25  k2  =  k3  —  /c,  j  —  2  (80) 

for  which  the  only  solution  is  k  =  1,  5  =  0.  The  boundary  condition  is  unstable  to  pertur¬ 
bations  away  from  the  unit  circle  and,  therefore,  exhibits  a  generalized  eigenvalue  instability 
at  the  boundary.  It  is  apparent  that  asymptotic  stability  for  the  semi-discrete  problem  does 
not  guarantee  Lax  stability. 

In  addition  to  the  previously  mentioned  examples,  Reddy  and  Treiethen  [14]  have  shown 
that  it  is  not  sufficient  to  consider  the  exact  eigenvalues  in  determining  the  stability  of  a 
method.  The  famous  Kreiss  matrix  theorem  gives  necessary  and  sufficient  conditions  for 
Lax  stability  in  terms  of  the  eigenvalues  of  the  matrix  M.  A  useful  and  equivalent  test  for 
determining  stability  is  the  analysis  of  the  resolvent  condition,  which  Reddy  and  Trefethen 
interpreted  as  involving  not  only  the  eigenvalues  of  the  matrix  M  but  also  the  e-pseudo- 
spectrum  of  the  discretization  matrix.  This  pseudo-spectrum  is  obtained  by  perturbing 
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the  matrix  M  by  an  arbitrary  matrix  of  norm  t.  Examples  where  the  e-pseudo-eigenvalues 
determine  the  stability  bounds  of  numerical  methods  are  given  by  Reddy  et  al.  [14].  It  is 
apparent  that  basing  a  stability  definition  on  the  eigenvalue  determination  of  the  spatial 
operator  is  not  sufficient  for  stability.  It  is,  however,  useful  for  the  present  applications  of 
higher-order  scheme  to  restrict  the  allowable  numerical  discretizations  to  those  possessing 
Lax  stability  and  the  property  of  bounded  LII-P  eigenvalues.  For  a  broad  class  of  spatial 
discretizations,  these  constraints  are  sufficient  for  the  stability  of  the  resulting  numeri  al 
scheme.  It  is  concluded  that  constraining  the  spatial  operators  to  those  possessing  these 
properties  is  a  useful  and  fairly  general  enhancement,  and  will  be  pursued  in  the  remainder 
of  this  work. 

Before  leaving  fourth-order  spatial  discretizations,  it  is  desirable  to  devise  a  uniformly 
fourth-order  scheme.  We  already  know  that  conventional  discretization  formulas  at  the 
boundaries  result  in  G-K-S  stable,  but  not  asymptotically  stable  schemes  for  both  the  explicit 
and  compact  cases.  The  NBS’s  used  in  each  case  relied  on  optimal  order  schemes  at  the 
boundary,  where  N+l  constraints  were  used  to  devise  the  Nth  order  scheme.  Specifically, 
five  grid-points  in  the  explicit  case,  and  four  grid-points  and  one  derivative  condition  in  the 
compact  case.  If  one  removes  the  constraint  of  using  optimal  schemes  at  the  boundaries, 
NBS  with  different  dissipative  characteristics  can  be  devised  and  an  asymptotically  stable 
spatial  scheme  which  is  uniformly  fourth-order  can  be  found. 

We  begin  by  devising  an  asymptotically  stable  fourth-order  compact  scheme,  which  we 
shall  denote  as  (43-4-43).  The  “43”  signifies  that  the  boundary  point  is  closed  with  a  fourth- 
order  three-parameter  family  of  schemes.  The  scheme  defined  at  grid-point  j=0  can  be 
written  as 


-7T -  =  (CoUo  +  CJA  +  C2U2  +  C3U3  +  C4U4  +  CdU  +  CeUe  +  C7U7)/(Ax).  (81) 

ax 

To  be  formally  fourth-order  accurate,  Taylor  series  truncation  analysis  relates  the  coefficients 
in  the  following  manner 

Co  =  -(a  -  28/4  +  3227  +  13068)/5040  (82) 

C'l  =  +(q  -  27/4  +  2957  +  5040)/720 

C2  =  — (q  —  26/4  +  2707  +  2520 ) /240 

C3  =  +(o  -  25:4  +  2477  +  1680)/ 144 

64  =  -(a  -  24b  +  2267  +  1260) / 144 

C5  =  +  (a  -  23/b  +  2077  +  1008)/ 240 

C6  =  -(a  -  22b  +  1907  +  840)/720 

C7  =  +  (o  -  21b  +  1 757  +  720) /504 0 


with  similar  expressions  defined  for  the  closure  at  the  other  end  of  the  domain.  By  sys¬ 
tematically  searching  the  three-parameter  space  spanned  by  the  parameters  a,/?,  and  7, 
coefficients  can  be  obtained  which  yield  an  eigenvalue  spectrum  which  is  bounded  to  the 
LH-P  The  values  of  a ,  /if ,  and  7  are  not  unique,  and  no  attempt  to  optimize  the  spectrum 
was  made.  A  particular  set  of  coefficients  which  make  the  scheme  asymptotically  stable  are 
a  =  — 1560,  3  =  —355  and  7  =  —35.  Figure  (9)  shows  the  resulting  spectrum  from  the 
(2-4-2),  (3-4-3)  and  (43-4-43)  schemes.  In  all  cases  the  eigenvalues  are  bounded  to  the  LH-P, 
and  the  resulting  scheme  is  asymptotically  stable. 

Since  the  asymptotic  stability  condition  'R.(Sj)  <  u  for  u  =  0  is  a  very  strong  necessary 
condition,  but  not  a  sufficient  condition,  we  still  must  show  G-K-S  stability  for  this  case. 
For  the  outflow  problem  the  model  equation  is 


dU_ 

dt 


dU_ 

dx 


=  0,  a-  >  0,  t  >  0; 


(83) 


no  boundary  condition  is  required  in  this  problem,  although  a  NBS  is  imposed  at  x  =  0. 
Assuming  a  solution  of  the  form  Vj(t)  =  exp5*  <?j,  where  <pj  —  <J>qk: ,  and  substitution  into 
the  inner  scheme  produces  the  resolvent  condition  for  the  eigenvalue  S 

( — f  4  +  k)5  =  3(  h  k)  (84) 

K  K 

at  each  grid-point  j  >  1.  At  grid-point  j=0,  the  parameter  scheme  produces  an  equation  of 
the  form 


3605  =  —(+727  —  1370k  +  1110k2  —  875k3  -F  775k4  —  552k5  +  220k6  —  35/c7).  (85) 

Solving  the  resolvent  condition  for  S  and  substituting  into  the  boundary  scheme  yields  a 
polynomial  in  k  of  the  form 

(k  -  1)5(35k4  +  95k3  -  168k2  -  227 k  -  353)  =  0.  (86) 


Solving  the  polynomial  for  k  produces  no  roots  which  are  in  magnitude  less  than  one:  thus, 
a  stable  condition.  The  possibility  of  k  =  1  exists  and  must  be  checked  for  generalized 
eigenvalues.  The  condition  is  the  same  one  tested  for  outflow  stability  in  eqn  (51),  and  was 
shown  to  be  stable.  Thus,  the  parametric  fourth-order  outflow  scheme  is  G-K-S  stable  for 
the  parameters  a.  ii  and  7  presented  above. 

To  show  stability  of  the  parameter  scheme  at  the  inflow,  we  study  the  partial  differential 
equation 


dV 

~dt 


+  _  =  0,  x  >  0,  /  >  0; 

ax 


(87) 
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with  the  boundary  condition  imposed  at  x  =  0  of  the  form 


U(0,t)=g(t),  t  >  0. 


(88) 


Eliminating  ^  between  the  boundary  scheme  at  j=0  and  the  inner-scheme  at  j  =  l,  yields  a 
combination  boundary  scheme  of  the  form 


1410 


dUx 

dx 


+  360 


dU2 

dx 


-(+353 74  +  137074  -  2190 U2  +  875 74  -  77574 

(89) 


+55274  -  22074  +  3574)/(Ax) 


Substituting  \j(t)  =  exp St  <j>0KJ  into  the  inner-scheme  and  the  boundary  scheme  produces 
the  resolvent  condition  eqn  (42),  and 


(1440k  +  360k2)5 


(+35374  +  1370k1  -  2190k2  +  875k3  -  775k4  +  552k5 


-220k6  +  35k7). 


(90) 


Solving  eqn  (42)  for  S  and  substituting  into  eqn  (91),  with  the  condition  that  7/0  —  0  yields 
a  polynomial  in  k  of  the  form 


k(— 2950  +  2210k1  -  2195k2  +  1615k3  -  1673k4  +  1213k5  -  293k6  -  80k7  +  35k8)  =  0.  (91) 


Solving  this  polynomial  results  in  no  roots  for  which  |k|  <  1 .  The  only  possibility  for  instabil¬ 
ity  is  the  condition  |k|  =  1.  Again,  this  condition  was  checked  previously  for  the  other  fourth- 
order  schemes,  and  was  shown  to  be  stable  for  the  inflow.  The  compact  three-parameter 
family  (43-4-43)  is  stable  for  inflow  and  outflow,  if  the  parameters  are  a  —  —1560,/?  =  —355 
and  7  =  —35.  In  addition,  the  scheme  produces  an  eigenvalue  spectrum  for  the  scalar  wave 
equation  which  is  bounded  to  the  Left  Half-Plane.  It  can  be  shown  that  treating  the  inflow 
and  outflow  boundaries  of  the  explicit  fourth-order  scheme  in  a  similar  manner  (43, 4-4-4, 43) 
also  produces  an  asymptotically  stable  scheme.  Whether  the  resulting  scheme  is  G-K-S 
stable  was  not  checked  in  this  work. 


8.  SIXTH-ORDER  SCHEMES 


As  a  last  step  in  this  work,  the  ideas  and  techniques  used  to  analyze  the  fourth-order 
compact  schemes  are  applied  to  sixth-order  compact  schemes.  All  the  schemes  tested  here 
will  be  based  on  the  sixth-order  compact  inner-scheme  developed  by  Lele  [15].  The  scheme 
can  be  written  in  the  form 

—  74-2  —  287  4— i  T  2874+1  +  74+2 


<97/,-i  +  3<974  +  7774+, 


dx 


dx 


dx 


12  Ax 


9=2,  ... ,  TV  —  2.  (92) 
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The  scheme  utilizes  information  from  five-point  explicitly  and  three  point  implicitly.  As  a 
consequence  of  the  five  point  width  of  the  stencil,  NBS’s  must  be  provided  at  two  points  at 
each  end  of  the  domain  (j=0,l  and  j=N-l,N).  The  physical  boundary  condition  is  used  at 
one  of  the  inflow  points.  To  ensure  formal  sixth-order  accuracy  for  the  hyperbolic  problem, 
the  boundary  points  must  be  closed  with  at  least  fifth-order  formulas,  which  for  the  optimal 
schemes  the  short  hand  nomenclature  would  be  (5, 5-6-5, 5).  (Again,  note  that  there  is  no 
ambiguity  in  comparing  this  nomenclature  with  that  from  the  sixth-order  explicit  scheme, 
since  it  would  require  three  boundary  formulas  at  each  end  of  the  domain.)  In  keeping  with 
the  convention  of  this  work,  the  closure  at  each  end  of  the  domain  is  done  in  a  asymmetric 
manner,  so  that  either  the  inflow  or  the  outflow  problem  can  be  easily  accommodated. 

It  has  proven  extremely  difficult  to  find  G-K-S  stable  schemes  which  are  formally  sixth- 
order  accurate.  We,  therefore,  begin  the  discussion  by  presenting  the  stability  analysis 
of  a  family  of  lower-order  schemes.  Two  of  the  schemes  in  this  family  are  the  1)  (3,5-6- 
5,3)  and  2)  (4, 5-6-5, 4).  The  formal  accuracies  of  these  schemes  are  fourth-  and  fifth-order, 
respectively.  The  closure  at  the  grid-point  j=0  (with  corresponding  formulas  written  at  j=N) 
is  accomplished  by 

— oUo  +  4(/i  +  l(2 


dU0  2dU  1 


dx 

dUo 

dx 


+  3 


dx 

dU\ 

dx 


2  Ax 

—  17 (Jo  +  9U\  4-  9U2  —  U3 
6Ax 


(93) 

(94) 


for  the  third-  and  fourth-order  schemes,  respectively,  while  the  fifth-order  closure  at  the 
point  next  to  the  wall  is  accomplished  in  all  cases  by  the  scheme 

—  lOf'o  —  9  U\  4-  I8U2  +  U3 


dU0  dU 1  dU2 

+  6-77—  +  3- 


(95) 


dx  dx  dx  3A.r 

Wider  spatial  stencils  produce  stability  polynomials  of  dramatically  increased  complexity. 
Despite  the  fact  that  MACSYMA  was  used  to  determine  all  of  the  spatial  formulas  and 
the  stability  polynomial  of  the  sixth-order  schemes,  the  possibility  for  error  still  exists. 
Noting  the  ability  of  the  G-K-S  theory  to  accurately  predict  the  stability  envelop  of  the 
one- parameter  family  of  fourth-order  compact  schemes  ( 1 1  -4-4 ) ,  a  simple  test  was  devised 
to  verify  the  accuracy  of  the  G-K-S  calculations.  A  one-parameter  family  of  schemes  was 
created  by  combining  the  third-  and  fourth-order  closure  formula  at  each  end  of  the  domain. 
Symbolically,  the  combined  scheme  is  represented  by  (3l, 5-6-5, 31),  and  is  written  as 


jdlJ 0  ,^Ui  —  7)Uo  +  4f 1  +  T2l  \\dUo  ,  r,9Ui 

nl17  +  2~0z - 251 - 1  + 11  “  a)l  a7  + 3  aT 


(96) 


1 7 f'o  +  9b  1  4-  9U2  —  U3 
6Ax 


=  0. 
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For  a  =  0  (or  1),  the  scheme  produces  the  optimal  fourth-order  (or  third-order)  variant,  and 
for  all  other  values  of  a,  the  formula  is  a  third-order  scheme. 

The  model  equation  for  the  outflow  quarter-plane  problem  is  the  same  as  described  in 
eqn  (83).  Solutions  of  the  form  Uj(t)  =  exp5t  cf)0KJ ,  satisfy  the  numerical  scheme,  giving  the 
sixth-order  inner-scheme  the  resolvent  condition 

1  1  oo 

( — F  3  +  k)S  =  ( - 2 - f  28/e  +  ft2)/ 12  (97) 

K  K  K 

for  the  eigenvalue  S.  There  are,  in  general,  two  roots  to  the  resolvent  equation  which  will 
yield  |k|  <  1  for  large  'R(S),  since  the  resolvent  is  a  polynomial  of  degree  four.  The  other 
two  roots  will  become  exponentially  unbounded  as  j  becomes  large  and  can  be  ignored.  The 
general  solution  has  the  form  Uj(t)  =  exp^Cj/cF  +  C2kF).  Substituting  this  expression  into 
the  boundary  eqns  (95,97)  yields  expressions  for  the  constants  Cx  and  C2.  The  expressions 
are 


C'iT’o(fti)  +  C2Fq(k2)  —  0 

Ci-Fi(fti)  +  C2Fi(k2)  =  0  (98) 


where 

F0(ft)  =  ((  —  (6c*  +  18)/c  -  18ft2)S  -  ((2a  +  3)  +  (3a  +  27)«  -  (6a  +  27)k2  +  (a  -  3)ft3)) 
Fi(«)  =  ((1  +  6k  +  3«2)S-  (-10 -9ft  +  18ft2  +  lft3)/3).  (99) 

Equation  (98)  has  only  the  trivial  solution  unless  the  determinant  condition  F0(fti)Fi(«:2)  — 
Fi(ft2)F0(fti)  —  0  is  satisfied.  Solving  the  resolvent  eqn  (97)  for  S  and  substituting  into  the 
determinant  condition  yields  an  expression  relating  the  two  «’s  as 

(«j  —  l)3(ft2  —  l)3((7a  —  3)«ik2  +  (—2a  +  3)(«i  +  k2)  —  3a  —  3)  =  0.  (100) 


Solving  eqn  (100)  yields  Kq  =  1  or  n2  =  1,  or 


(2a  —  3)k2  +  3a  +  3 
(7a  —  3)«2  —  2a  +  3 


(101) 


To  obtain  an  additional  independent  relationship  for  ft]  and  k2,  the  resolvent  condition  eqn 
(97)  is  used.  Solving  the  resolvent  equation  for  S  and  noting  that  both  K\  and  k2  satisfy 
this  expression  for  S,  yields 


-  -7  +  28ft i  +  ft/2) 


S’  = 


*1 


~  +  28k2  +  K22) 


12(^7  +  3  +  fti) 


12(^-  +  3  +  k2) 


(102) 
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Combining  eqn  (101)  and  eqn  (102)  produces  a  single  sixth-order  polynomial  in  x2,  for  which 
the  roots  can  be  found  numerically.  Equation  (101)  and  eqn  (97)  are  then  used  to  determine 
the  numerical  values  of  xt  and  S.  An  eigensolution  exists  for  the  problem  if,  for  |/cj|  <  1  and 
|x2|  <  1,  there  exists  a  S  with  real  part  greater  than  zero.  Solving  the  outflow  polynomial 
for  x2  yields  the  result  that  the  boundary  is  unstable  for  —9.16  <  a  <  —1.86.  The  two 
limiting  cases  (a  =  0  and  a  =  1)  for  which  the  scheme  is  fourth-  or  third-order  accurate,  are 
stable.  It  should  be  noted  that  the  possibility  of  instability  also  exists  for  the  case  —  x2, 
and  for  xq  or  x2  =  1.  None  of  these  conditions  showed  instability,  however. 

The  model  equation  for  the  inflow  quarter-plane  problem  is  the  same  as  described  in  eqns 
(87,88).  Solutions  of  the  form  Uj(t)  —  expSt  4>qkj ,  will  satisfy  the  numerical  scheme,  giving 
the  sixth-order  inner-scheme  the  resolvent  condition 

( — b  3  +  x)S  —  — ( — —  — - b  28 x  +  x2)/ 12.  ( 1 03 ) 

X  Kz  X 

Eliminating  between  the  boundary  schemes  at  grid-points  j=0  and  j  =  l,  yields  an  ex¬ 
pression  written  for  grid-point  j  =  l  of  the  form 

(6a  +  18)^- +  18^  =  ((2a  +  3)£/0  +  (3a  +  27)£/i -(6a  +  27)tfi 
Ox  ox 

(104) 


+(a  -  3)t/3))/(6Ax). 


The  general  solution  has  the  form  Uj(t)  =  exp5t(6'ix1-'  +  C2x2ji).  Substituting  this  expression 
into  the  boundary  expressions  at  j  =  1 ,2  yields  two  expressions  for  the  constants  C\  and  C2. 
The  expressions  are 


C1F1(x1)  +  C2F,(x2)  =  0 

C1F2(x1)  +  C2F2(x2)  =  0  (105) 


where 

Fi(k)  =  ((-(6a  +  18)x  -  18x2)5  -  (+(3a  +  27)x  -  (6a  +  27)x2  +  (a  -  3)x3)) 

F2(x)  =  1.  (106) 

The  simple  form  of  the  expression  F2  =  1  results  from  reducing  the  modal  equation  at  point 
j=2  (with  U0  set  to  zero)  by  use  of  the  resolvent  condition  (with  Uq  not  equal  to  zero). 
Equation  (105)  can  only  have  a  nontrivial  solution  if  the  determinant  is  identically  zero. 
Solving  eqn  (103)  for  5  and  substituting  into  the  determinant  condition  from  eqn  (105) 
yields 

(2a  -  3)xt5  +  (-5a  +  15)x1‘l  -  30x-i3  +  (6a  +  24)x,2  +  (-22a  -  33)x11  -  (a  +  3) 

2x-i2  +  6xq 1  -f  2 
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(2a  -  3)k25  +  (-5a  +  15)k24  -  3(k23  +  (6a  +  24)/c22  +  (-22a  -  33)k21  -  (a  +  3) 

2k22  +  +  2 


0. 


(107) 

Together  with  eqn  (102),  they  provide  two  equations  for  the  unknown  «q  and  n2. 

The  inflow  polynomial  equations  are  far  more  difficult  to  solve  since  they  do  not  factor 
appreciably.  A  change  of  variables  from  K\  and  k2  to  “x”  and  “y”,  simplifies  the  algebra. 
Substituting 


T  ^2  —  2 y 

K\K2  =  x 


into  eqn  (102)  yields 

(— 8x  -  8 )y2  +  (  — 12x2  -  224x  -  12 )y  -  2x3  -  166x2  -  166x  -  2  =  0 


(108) 


(109) 


which  has  the  solution 

—  (3x3  -f  56x  +  3  ±  \/5x4  +  2490x2  +  5) 
y  =  - - 

for  x  7^  —1.  (The  case  x  =  -1  degenerates  into  y=0,  for  which  «q  =  —  k2  =  ±1,  a  condition 
which  produces  no  eigensolutions).  Either  root  can  be  used  since  the  final  polynomial  results 
from  squaring  an  intermediate  result  to  clear  the  square  root  in  the  expression.  Substituting 
eqn  (108)  into  eqn  (107)  and  further  simplifying  with  the  expression  for  y  from  eqn  (110) 
yields  a  ninth-order  polynomial  in  the  variable  x,  with  coefficients  that  are  functions  of  the 
variable  a.  Solving  this  expression  numerically  yields  the  roots  for  x.  The  values  of  y  are 
then  determined  from  eqn  (110)  and  the  values  of  «q  and  k2  are  obtained  from  eqn  (108). 
Again,  the  condition  k  =  1  and  the  case  Kq  =  k2  did  not  show  instability  over  the  parameter 
range  tested  in  this  study.  For  the  inflow  boundary  condition,  the  instability  envelop  for  the 
parameter  a  was  determined  to  be  —1.86  <  a  <  —0.447.  It  is  evident  that  the  two  degenerate 
conditions,  a  =  0  and  a  =  1,  are  G-K-S  stable  for  the  inflow.  By  combining  the  inflow  and 
outflow  results,  the  theoretical  (G-K-S)  range  of  instability  for  the  one-parameter  family  of 
schemes  (31 ,5-6-5, 31 ),  was  determined  to  be  —9.18  <  a  <  —0.447.  The  two  degenerate  cases, 
(3, 5-6-5, 3)  and  (4, 5-6-5, 4),  were  both  G-K-S  stable  schemes. 

To  determine  the  accuracy  of  the  G-K-S  calculations,  another  method  of  showing  stability 
was  used,  arid  the  results  were  compared.  The  necessary  condition  for  stability  on  the 
eigenvalue  structure  7 Z(Smax)  provides  such  a  test.  Figure  (10)  shows  the  results  of  an 
eigenvalue  determination  spanning  the  range  — 10  <  a  <  2,  as  determined  numerically.  The 
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maximum  eigenvalue  7Z(Smax )  (the  eigenvalue  with  the  largest  absolute  real  component), 
is  plotted  as  a  function  of  the  parameter  a  for  three  grid  densities.  Since  Lax  stability  in 
a  finite  domain  requires  that  7 Z(S)  <  u>  for  u>  >  0,  the  7 Z(Smax)  should  remain  bounded 
with  increasing  grid  density  if  the  scheme  is  to  be  stable  for  a  particular  value  of  a.  It  is 
apparent  from  Figure  (10),  that  the  stability  boundary  at  a  =  —9.15  is  accurately  predicted 
by  G-K-S  theory.  The  stability  boundary  at  a  —  —0.45  is  less  well  defined  in  the  eigenvalue 
determination,  and  must  be  further  investigated  to  show  the  correlation  between  G-K-S 
theory  and  eigenvalue  determination.  On  the  relatively  coarse  grids  presented  in  Figure 
(10),  the  maximum  eigenvalues  near  the  limit  a  =  —0.45  are  growing  with  increasing  grid 
density.  Whether  they  are  growing  in  a  bounded  manner  determines  if  they  satisfy  the 
necessary  condition  for  stability.  Table  (2)  shows  the  behavior  of  TZ{Smax)  for  various  grid 
densities  at  a  =  —0.40. 


Grid 

1Z{Smax) 

A% 

21 

1.000 

NA 

41 

1.250 

25 

81 

1.360 

8 

161 

1.582 

16 

321 

1.704 

8 

641 

1.780 

5 

Table  2:  Grid  convergence  of  7 Z(Smax)  for  a  =  —0.40. 

As  the  number  of  grid-points  becomes  larger,  the  maximum  eigenvalue  asymptotes  to  a 
constant  u>,  thus  a  stable  condition.  Note  that  the  convergence  to  the  asymptotic  limit 
is  slow  for  q  =  —0.40.  Similar  grid  refinement  studies  at  values  of  a  =  -0.45  and  -0.50 
showed  linear  growth  for  all  grids,  thus  an  unstable  condition.  Based  on  a  numerically 
determ. ned  eigenvalue  determination  over  the  range  of  —10  <  a  <  2,  the  G-K-S  theory  is 
accurately  predicting  the  stability  envelop.  It  should  be  noted  that  the  slow  convergence 
to  the  asymptotic  limit  (and  the  fact  that  the  test  is  a  necessary  not  sufficient  one)  makes 
testing  for  stability  by  numerical  eigenvalue  determination  somewhat  unreliable.  For  many 
non-borderline  cases,  it  does  provide  an  accurate  measure  of  stability. 

The  eigenvalue  determination  provides  information  on  the  asymptotic  stability  of  the 
schemes  as  well.  If  for  all  the  grids,  the  TZ(Smax)  <  0,  the  scheme  is  asymptotically  stable. 
For  valves  of  the  parameter  a  <  —9.15  and  a  >  0.4,  eigenvalue  determination  indicates 
asymptotic  as  well  as  Lax  stability  of  the  resulting  scheme.  Figure  (11)  shows  the  eigenvalue 
spectrum  of  the  (3, 5-6-5, 3)  and  (4, 5-6-5, 4)  schemes.  It  is  apparent  that  both  satisfy  the 
necessary  condition  of  Lax  stability,  and  that  the  (3, 5-6-5, 3)  scheme  is  asymptotically  stable. 
Because  of  the  eigenvalues  in  the  R II- P,  the  (4, 5-6-5, 4)  scheme  does  not  exhibit  asymptotic 


stability.  Figure  (12)  shows  a  plot  of  the  error  of  the  solution  to  the  scalar  wave  equation 
defined  by  eqns  (28,  29,  30)  when  discretized  with  the  (4, 5-6-5, 4)  scheme.  The  time  interval 
was  0  <  t  <  100  ,  and  time  was  advanced  with  a  fourth-order  Runge-Kutta  scheme.  The 
logl0  of  the  Li  error  is  plotted  as  a  function  of  time  for  three  grid  densities:  21,  41,  and  81 
points,  respectively.  For  a  CFL  =  1,  the  error  does  not  grow  in  time.  For  CFL’s  of  order 
0.1,  the  error  is  nearly  uniform  for  a  period  of  time,  then  begins  to  grow  exponentially  with 
time.  In  all  cases,  doubling  the  grid  decreases  the  error  of  the  simulation  by  a  factor  of  16 
-  32  (error  is  dominated  by  either  the  fourth-order  time,  or  the  fifth-order  space  truncation 
terms.)  The  amplification  is  accurately  predicted  by  the  eigenvalue  determination.  Table  (3) 
shows  the  numerical  amplification  rate  compared  with  the  /R-(Smax)  ,  for  which  the  agreement 
is  excellent. 


Grid 

Numerical 

QTC(5mal) 

21 

0.1245 

0.1228 

41 

0.1402 

0.1381 

81 

0.1351 

0.1354 

Table  3:  Numerical  vs.  Theoretical  Growth  Rate;  (4, 5-6-5. 4). 

The  preceding  examples  of  sixth-order  schemes  show  the  strength  of  G-K-S  analysis  to  ac¬ 
curately  predict  the  stability  of  complex  higher-order  schemes.  They  also  show  the  intimate 
relationship  between  the  eigenvalues  of  the  spatial  operator  and  the  stability  of  the  resulting 
scheme. 

We  now  present  schemes  that  are  formally  sixth-order  accurate,  being  closed  at  the 
boundaries  by  at  least  fifth-order  stencils.  Our  first  attempt  is  with  optimal  fifth-order 
closure  at  the  boundaries,  resulting  in  the  scheme  (5, 5-6-5, 5).  Figure  (13)  shows  the  error 
of  the  simulation  of  the  scalar  wave  equation  defined  by  eqns  (28,  29,  30).  The  behavior  of 
this  scheme  is  fundamentally  different  from  the  (4, 5-6-5, 4)  scheme  in  several  ways.  On  all 
grids,  the  error  is  always  monotonically  increasing  in  time.  For  CFL’s  near  the  theoretical 
maximum  value,  the  error  increases  at  a  lower  rate,  but  is  not  suppressed  as  it  was  with 
the  lower-order  schemes.  In  addition,  the  exponential  growth  rate  of  the  error  increases 
with  increasing  grid  density.  On  the  grids  shown,  the  error  in  the  solution  could  not  be 
systematically  reduced  by  refining  the  grid  and  repeating  the  calculation  to  a  specified  time 
level  T.  In  spite  of  these  differences,  the  eigenvalue  determination  still  accurately  predicts 
the  growth  of  the  solution.  Table  (4)  shows  a  comparison  of  the  numerical  and  theoretical 
amplification  rates.  The  theoretical  values  are  determined  from  the  TZ(Smnr)  f°r  each  grid. 
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Grid 

^ Numerical 

Q'7?(S„laj) 

21 

0.7121 

0.7138 

41 

1.104 

1.010 

81 

1.749 

1.742 

Table  4:  Numerical  vs.  Theoretical  Growth  Rate;  (5, 5-6-5, 5). 

Again,  the  agreement  is  excellent.  The  solution  grows  at  a  rate  which  for  long  times  is 
dominated  by  the  eigenvalue  with  the  maximum  real  part.  The  eigenvalue  determination 
accurately  predicts  the  behavior  of  the  numerical  solution  even  for  this  case  which  appears 
not  to  be  Lax  stable. 

Figure  (14)  shows  a  plot  of  the  eigenvalue  spectrum  for  the  (5, 5-6-5, 5)  scheme  on  a 
21,  41  and  81  grid.  The  7?.(S'mar)  is  obviously  increasing  for  these  grids,  and  appears  to 
be  increasing  without  bound,  as  opposed  to  an  asymptotic  limit.  This  wouH  violate  a 
necessary  condition  for  Lax  stability.  As  was  seen  in  the  (31, 5-6-5, 31)  example,  one  cannot 
draw  a  precise  conclusion  from  grid  refined  eigenvalue  determination,  although,  the  trends 
have  in  the  cases  presented  thus  far  been  the  same.  We  must  ultimately  rely  on  G-K-S 
stability  theory  for  the  final  answer  as  to  whether  the  scheme  is  stable. 

We  begin  by  determining  the  stability  cf  the  outflow  boundary  for  the  (5, 5-6-5, 5)  scheme. 
The  quarter-plane  problem  appropriate  for  this  analysis  is  that  described  in  eqn  (83).  No 
boundary  conditions  are  necessary,  but  NBS’s  are  used  at  grid-points  j =0,1 .  The  closure  at 
grid-point  0  is  accomplished  with 

dUo  dU i  — 37f70  -f  SUi  +  366-^  —  8U3  +  U4 

-o7  +  iH7  =  - I2S -  (m) 

while  that  of  grid-point  j  =  l  is  accomplished  with  eqn  (95).  Solutions  of  the  form  Uj(t)  = 
exp6'  0q will  satisfy  the  numerical  scheme,  giving  the  sixth-order  inner-scheme  the  resol¬ 
vent  condition  shown  in  eqn  (97).  1  here  will  be  two  roots  to  the  fourth-order  polynomial  in  k 
which  are  |k|  <  1,  and  the  general  solution  will  have  the  form  Uj(t)  =  exp5t(Ci«:1J  +  C^^2)- 
Substituting  this  expression  into  the  two  boundary  conditions  yields  boundary  expressions 
for  the  constants  C\  and  C 2  of  the  form 


CiFo(«i)  +  C-2  Fo(k2)  —  0 

Ci/*i(aci)  +  C^Fi^)  =  0  (112) 


where 

F0(k)  =  ((1 +4^)5- (-37  +  8k  +  36k2 -Sk3  1  k4)/12) 

Fj(k)  =  ((1  +6k  +  3k2)S’-(-10-9k  +  18k2  +  k3)/3).  (113) 
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=  0. 


(114) 


Equation  (112)  can  only  have  a  nontrivial  solution  if  the  determinant  condition  Fo(kx  )F\(k2)  — 
F\(k2)Fo(i<i)  =  0  is  satisfied.  Solving  the  resolvent  condition  eqn(97)  for  S  and  substituting 
into  the  determinant  condition  produces 

(fv'i  —  1)4(k2  —  1)  {k2  ~  ki)(4k2ki  +  *2  +  /Cl  —  6)  _ 

144 /Cl  K2(k^2  — 3 AC  1  +  1  )(k22  +  3k2  4“  1) 

Solving  eqn  (114)  yields  Ki  =  1  or  «2  =  1,  or  kx  =  k2,  or  the  expression 

k,  =  -JH2Z±.  (115) 

4k2  +  1 

The  first  two  roots  are  the  same  roots  that  have  been  shown  previously  to  be  stable.  The 
condition  that  both  roots  must  be  |/c|  <  1  precludes  the  last  root.  The  third  root  can 
be  shown  to  be  stable  by  testing  the  derivative  condition  of  the  polynomial  as  follows. 
Multiplying  and  dividing  eqn  (112)  by  the  non-zero  rows  and  columns  does  not  change  the 
roots  of  the  determinant  conditions.  The  resulting  expression  is 

Eo(k2)  —  F0(kx) 


F0{ki ) 


Fi(k\) 


K2  ~  «i 

Fi{k2)  -  F-i(ki) 


=  0. 


K-  2  ~  K 1 

Taking  the  limit  as  kx  — +  k2  yields  the  expression  for  the  determinant  condition  Fo 
Fi{n)df°^  =  0  .  The  resulting  expression  for  k  is 

,12 

=  0.  (lit)) 


(«-!)' 


144k2(k:2  +  3k  +  1)‘ 

The  outflow  boundary  is,  thus,  stable  for  the  (5, 5-6-5, 5)  scheme. 

The  model  equation  for  the  inflow  quarter-plane  problem  is  the  same  as  described  in 
eqn  (87,  88).  Solutions  of  the  form  Uj(t)  =  exp St<t>oKj,  will  satisfy  the  numerical  scheme, 
giving  the  sixth-order  inner-scheme  the  same  resolvent  condition  as  was  given  in  eqn  (103). 
Again,  the  general  solution  will  have  the  form  Uj{t)  =  expS((CiKiJ  +  C2K22).  Substituting 
this  expression  into  the  two  boundary  conditions  and  making  the  simplification  that  Uo  =  0 
gives  two  equations  for  the  constants  C\  and  C2  of  the  form 


where 


C\F\{K\)  +  C'2Fi(k2)  —  0 

C’i/'2(«l)  +  E*2  F2[  K*2  )  =  0 

F\ ( k )  =  ((2k  +  3k2) 5  -  (  1  14k  -  36k2  -  12k3  +  k4)/1‘2 
F2(k)  =  1. 


(117) 


(118) 
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The  determinant  condition  for  which  there  is  a  nontrivial  solution  simplifies  to  Fl(/ci)  — 
F\{k2)  =  0  and,  along  with  eqn  (102),  provide  two  equations  for  the  two  unknowns  and 
k2.  Using  the  change  of  variables  aq  +  k2  =  2 y  and  k^k2  —  x,  eqn  (118)  becomes 

—64 y5  +  (192  -  96 x)y4  +  (-16x2  +  352a.-  -  240)*/3  +  (120x2  -  504x  +  160)y2 

+  (8:r3  -  216a:2  +  360x  -  56)y  -  18x3  +  142x2  -  142x  +  18  =  0.  (119) 

Substituting  the  functional  relationship  y=y(x)  provided  by  eqn  (110)  into  eqn  (119)  and 
simplifying  yields 

(x  -  1 ) (x  +  l)5(x12  +  261xn  +  24298x10  +  864903x9  +  864903x9  +  5558711x8 
+  1650241  Ox7  +  27479264x6  +  2853S822x5  +  5255107x4  +  429169x3 


+  18614x2  +  435x  +  5)  =  0. 


(120) 


Two  mots  to  this  polynomial  give  rise  to  eigenvalues  S  which  are  in  the  RH-P.  They  are  x  = 

(-0.01969168445  ±0.023980223190,  which  yields  Kl  =  (0.157055621  +  0.9436011290,  with 
k2  =  (-0.02810826491, +0.016190234380  and  5  -  (0.0428389  ±  1.399440-  The  numerical 

solutions  satisfy  the  governing  equations  to  approximately  machine  precision  (1.0e-13).  Thus, 
the  inflow  for  the  (5, 5-6-5, 5)  scheme  is  G-K-S  unstable.  This  verifies  the  trends  presented 
earlier  in  the  eigenvalue  grid  refinement  analysis,  and  the  simulation  of  the  scalar  wave 
equation,  both  of  which  indicated  the  sixth-order  scheme  was  Lax  unstable. 

The  unstable  (5, 5-6-5, 5)  scheme  previously  discussed  was  implemented  using  optimal 
fifth-order  boundary  formulas.  By  relaxing  the  constraint  of  optimal  order  schemes  at  the 
boundaries,  the  possibility  of  a  fifth-order  closure  which  will  be  G-K-S  stable  still  exists. 
Taylor  series  truncation  analysis  was  used  to  develop  parametric  relations  for  the  closure 
formulas  at  the  two  NBS’s.  They  were  constrained  to  be  fifth-order,  and  explicit  in  nature. 
To  facilitate  a  wide  range  of  closures,  each  point  was  given  two  degrees  of  freedom.  The 
symbolic  formula  for  the  new  scheme  is,  thus,  (52,52-6-52,52).  The  scheme  defined  at  grid- 
point  j=0  and  1  can  be  written  as 


^  -  (C0 bo  ±  Cdh  +  C2u2  +  C3U3  +  C4U4  +  cw,  ±  C6u6  +  C7U7)/( Ax) 

ox 

— —  —  (  DqL  u  +  /2|f  1  +  D2U2  +  D,\ b- 3  +  1  +  D$V 5  +  b)6b6  +  D7U7) /(i\x) 

ox 


(121) 

(122) 


where 


Co  =  -  ( o0  -  28/fo  +  1 3068 )/504  0 


(123) 
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Cj  =  +(q0  -  27/?0  +  5040)/720 
C2  =  -(q0  -  26/?0  +  2520)/240 
C3  =  +(a0  —  25/?o  +  1680)  / 144 
C4  =  — (a0  —  24/?o  +  1260)/ 144 
C5  =  -f  (u0  ~  23/?o  +  1008)/240 
C6  =  — (q0  —  22/?0  +  840)/720 
C7  =  +(a0  -  21/?0  +  720)/5040 

and 

D0  =  — (c*i  -  21^  +  720) /5040  (124; 

A  =  +(ai  —20/?!  —  1 044) /720 
D2  =  -(a,  -  19A  -  720)/240 
Z)3  =  +(<*,  -  18/5,  -  360)/144 
D4  =  —(a,  —  17/3,  —  240)/144 
A  =  +(«i  -  16/?!  -  180)/240 
A  =  — (c*i  -  15/?!  -  144)/720 
Dr  =  +(a,  -  14/?,  -  120)/5040 

with  similar  expressions  defined  for  the  closure  at  the  other  end  of  the  domain.  By  systemat¬ 
ically  searching  the  four-parameter  space  spanned  by  the  parameters  a0,/?o,ai,  and  /?,  with 
an  eigenvalue  code,  an  arbitrary  eigenvalue  spectrum  can  be  approximated.  A  particular 
set  of  coefficients  for  which  the  eigenvalue  spectrum  is  bounded  to  the  Left  Half-Plane  is 
ao  =  1809.257, /?o  =  —65. 1944,  qi  =  —262.16  and  /?,  =  —26.6742.  The  values  are  not  unique 
and  no  attempt  has  been  made  to  find  optimal  values  of  these  coefficients.  The  eigenvalue 
spectrum  for  this  case  is  shown  on  Figure  (14).  Note  that  the  shape  of  the  spectrum  is  similar 
to  that  of  the  (5, 5-6-5, 5)  scheme,  but  that  the  TZ(Smax)  <  0  instead  of  increasing  without 
bound.  The  scheme  satisfies  the  necessary  condition  for  Lax  stability,  and  is  asymptotically 
stable  by  our  definitions. 

The  stability  analysis  for  the  two  quarter-plane  problems  involved  in  establishing  the 
G-K-S  stability  of  the  new  (52, 52-6-52, 52)  is  extremely  formidable.  It  pushes  MACSYMA  to 
the  limits  of  its  capabilities  on  present  machines.  In  addition,  128-bit  arithmetic  is  required 
to  ensure  precision  when  determining  the  roots  of  the  resulting  polynomials  in  “x”.  The 
scheme  has  been  shown  to  be  G-K-S  stable  for  the  parameters  given  above  for  the  inflow  and 
outflow  problems.  Thus  a  formally  sixth-order  scheme  has  been  developed  which  is  G-K-S 
(Lax)  stable,  and  asymptotically  stable  for  the  scalar  case. 
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r°  verify  its  accuracy,  Table  (5)  shows  a  grid  refinement  study  performed  with  the  new 
sixth-order  scheme.  The  model  problem  is  the  scalar  wave  equation  defined  by  eqns  (28,  29, 
30).  The  time  advancement  scheme  is  the  fourth-order  Runge-Kutta  algorithm,  with  a  CFL 
of  0.10.  Temporal  refinement  studies  were  performed  to  ensures  that  the  leading  error  terms 
on  all  grids  were  from  the  spatial  discretization  operator.  Listed  is  the  grid  and  its  error, 
and  the  slope  between  each  successive  refinement.  The  asymptotic  stability  of  the  spatial 
operator  ensure  that  the  solution  does  not  grow  exponentially  for  long  times. 


Grid 

L^error 

slope 

21 

-2.363 

NA 

31 

-3.582 

-6.9 

41 

-4.225 

-5.1 

51 

-4.724 

-5.1 

61 

-5.150 

-5.4 

81 

-5.849 

-5.6 

101 

-6.406 

-5.7 

121 

-6.867 

-5.8 

Table  5:  Grid  refinement  study  of  the  (52, 52-6-52, 52)  scheme. 

The  data  point  corresponding  to  21  grid-points  is  erroneous  because  the  grid  is  too  coarse 
for  the  scheme  to  exhibit  its  higher-order  properties.  Note  that  in  this  example  the  scheme 
becomes  at  least  fifth-order  accurate  at  approximately  10  grid-point  s/(27r  radians).  In  the 
limit  N  — *  oo  where  N  is  the  number  of  grid-points,  the  scheme  shows  a  slope  of  -6,  the 
formal  accuracy  of  the  inner-scheme.  Note  that  for  N  small,  the  four  points  which  are  treated 
with  fifth-order  accuracy  degrade  the  formal  accuracy  by  one  degree. 

An  asymptotically  stable  sixth-order  spatially  accurate  scheme  has  been  developed  for 
use  in  a  method-of-lines  discretization  of  a  hyperbolic  partial  differential  equation.  The 
eigenvalues  of  the  new  scheme  are  for  the  scalar  case  bounded  to  the  Left  Half-Plane  of  the 
complex  plane  for  all  N.  The  necessary  condition  for  Lax  stability  is,  therefore,  satisfied. 
In  addition,  the  scheme  has  been  shown  to  be  G-K-S  stable  for  the  combined  inflow  and 
outflow  quarter-plane  analysis,  and  is,  therefore,  formally  Lax  stable  for  the  scalar  case.  For 
the  hyperbolic  system  of  equations,  the  use  of  any  of  the  Lax  stable  schemes  presented  in  this 
work  guarantees  the  Lax  stability  of  the  resulting  spatial  discretization  if  the  boundaries  are 
imposed  in  characteristic  form.  The  concept  of  asymptotic  stability  does  not  carry  over  from 
the  scalar  case  to  the  system.  For  the  case  of  the  system  with  all  of  the  physical  eigenvalues  of 
the  same  sign,  the  asymptotic  stability  is  still  retained.  For  the  case  of  mixed  eigenvalues,  it 
is  easy  to  show  that  even  though  being  asymptotically  stable  for  the  scalar  case,  exponential 
growth  of  the  solution  may  occur  for  boundaries  that  are  imposed  in  characteristic  form. 
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Work  continues  in  this  area,  to  determine  a  stronger  necessary  condition  for  asymptotic 
stability  that  will  allow  the  use  of  scalar  analysis  to  determine  spatial  schemes  which  are 
asymptotically  stable  for  the  system. 

9.  CONCLUSIONS 

The  stability  characteristics  of  various  compact  fourth-  and  sixth-order  spatial  operators 
were  assessed  using  the  theory  of  Gustafsson,  Kreiss  and  Sundstrom  (G-K-S)  for  the  semi¬ 
discrete  IBVP.  The  class  of  central  difference  interior  schemes  with  asymmetrically  closed 
boundaries  was  analyzed.  Because  of  formal  accuracy  considerations,  those  schemes  with 
boundary  closures  of  at  least  (N  —  l)th  spatial  order  for  an  ( N)th  order  inner  scheme  were 
the  focus  of  the  work.  It  was  found  that  conventional  third-  or  fourth-order  boundary 
conditions,  when  coupled  with  the  fourth-order  compact  inner-scheme,  resulted  in  a  G-K-S 
stable  scheme.  For  the  sixth-order  compact  inner-scheme,  the  conventional  boundary  closures 
of  fifth-  and  higher-order  were  found  to  be  G-K-S  unstable.  Fourth-order  and  lower-order 
closure  formulas  were  found  to  be  G-K-S  stable.  These  results  were  then  generalized  to  the 
fully  discrete  case  using  a  recently  developed  theory  of  Kreiss,  which  states  that  under  weak 
constraints,  the  stability  of  the  semi-discrete  operator  implies  stability  of  the  fully  discrete 
operator  if  a  locally  stable  temporal  method  is  used. 

The  conventional  definition  of  stability  was  then  sharpened  to  include  only  those  spatial 
discretizations  that  are  asymptotically  stable  (bounded  Left  Half-Plane  eigenvalues).  Many 
of  the  higher-order  schemes  which  are  G-K-S  stable  were  found  to  not  be  asymptotically 
stable.  Fourth-order  boundary  conditions  were  found  to  be  asymptotically  unstable  for  the 
schemes,  tested,  specifically:  1)  (4-4-4),  and  2)  (4, 5-6-5, 4).  A  series  of  compact  fourth-  and 
sixth-order  schemes  which  were  both  asymptotically  and  G-K-S  stable  were  then  developed. 
The  constraint  of  optimal  accuracy  from  a  specific  number  of  constraints  was  abandoned, 
thus  enabling  several-parameter  family  boundary  closures  to  be  developed.  A  three  param¬ 
eter  family  uniformly  fourth-order  scheme  (43-4-43),  as  well  as  a  four-parameter  sixth-order 
scheme  with  fifth-order  boundaries  (52, 52-6-52, 52)  were  developed  which  were  asymptotically 
stable.  No  attempt  was  made  to  optimize  the  parameters. 

All  the  schemes  which  were  found  to  be  G-K-S  stable  were  subjected  to  extensive  com¬ 
parisons  between  the  G-K-S  stability  predictions,  semi-discrete  eigenvalue  determination  and 
numerical  simulations.  In  all  cases,  consistent  and  complementary  results  were  achieved  with 
all  three  methods.  In  addition,  it  was  shown  that  the  eigenvalue  determination  accurately 
predicted  the  exponential  divergence  of  the  solution  for  the  cases  which  were  not  asymptot¬ 
ically  stable.  Work  continues  in  developing  asymptotic  definitions  for  the  case  of  systems  of 
equations. 
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PADE  BOUNDARY  CONDITION  ACCURACY 


Fig-1  Comparative  study  of  the  effects  of  boundary  closure  on  the  formal  accuracy  of  the 
fourth-order  compact  inner  scheme. 
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Fig-2  Comparative  study  of  the  effects  of  outflow  boundary  closure  on  the  formal  accuracy 
of  the  fourth-order  compact  inner  scheme. 
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Fig-3  Comparative  study  of  the  effects  of  inflow  boundary  closure  on  the  formal  accuracy 
of  the  fourth-order  compact  inner  scheme. 
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Fig-4  Numerical  determination  of  the  L 2  error  at  a  time  level  “T”  resulting  from  a  one- 
parameter  family  of  closure  formulas  at  the  inflow  boundary.  The  inner  scheme  is  the 
fourth-order  Pade  scheme. 


EIGENVALUE  SPECTRUM  (4-th  LINEAR) 
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Fig-6  The  numerically  determined  eigenvalue  spectrum  resulting  from  the  explicit  fourth- 
order  inner  scheme,  closed  with  either  third-  or  fourth-order  boundary  schemes. 
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7  The  numerically  determined  eigenvalue  spectrum  resulting  from  the  fourth-order  Pade 
inner  scheme,  closed  with  eithe :  third-  or  fourth-order  boundary  schemes. 
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AMPLIFICATION  FACTOR  (4th  PADE  ) 


C-L 


Fig-8  The  maximum  numerical  amplification  factor  as  a  function  of  CFL,  plotted  for  the 
uniformly  fourth-order  Pade  scheme,  with  a  fourth-order  Runge-Kutta  time  advance¬ 
ment  scheme. 
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EIGENVALUE  SPECTRUM  (NTH  PADE) 
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Fig-9  The  numerically  determined  eigenvalue  spectrum  resulting  from  the  fourth-order  Pade 
inner  scheme,  closed  with  either  second-,  third-  or  newly  developed  fourth-order  bound¬ 
ary  schemes,  on  various  grids. 


50 


MAXIMUM  E!GFNVA:.Ut  (p, 5-6-5, p) 


alpha 


dg-10  The  maximum  real  part  of  the  numerically  determined  eigenvalue  spectrum,  result¬ 
ing  from  a  sixth-order  Pade  inner  scheme  and  a  one-parameter  family  of  third-order 
boundary  schemes.  The  ordinant  a  is  the  boundary  parameter. 
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EIGENVALUE  SDECTRUV  (6lh  PADE  i 
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Fig-1 1  The  numerically  determined  eigenvalue  spectrum  resulting  from  the  sixth-order  Pade 
inner  scheme,  closed  with  either  third-  or  fourth-order  boundary  schemes,  on  various 


ERROR  OF  PADE  SCHEME  (4,5-6-5.4) 


tig  1^  L he  L 2  error  as  a  function  of  time,  resulting  from  the  uniformly  sixth-order  Pade 
scheme  with  fourth-order  boundary  closure  (4, 5-6-5, 4),  for  various  CFL’s  and  grid 
densities. 
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ERROR  OF  PADE  SCHEME  (5, 5-6-5, 5) 


Fig— 13  The  L2  error  as  a  function  of  time,  resulting  from  the  uniformly  sixth-order  Pade 
scheme  with  fifth-order  boundary  closure  (5, 5-6-5, 5),  for  various  CFL’s  and  grid  den¬ 
sities. 
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EIGENVALUE  SPECTRUM  (6th  PADE) 
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Fig -14  T  he  numerically  determined  eigenvalue  spectrum  resulting  from  the  sixth-order  Pade 
inner  scheme,  closed  with  conventional  fifth-  or  newly  developed  fifth-order  boundary 
schemes,  on  various  grids. 
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